SPH Simulation of Interior and Exterior Flow Field Characteristics of Porous Media

: At the present time, one of the most relevant challenges in marine and ocean engineering and practice is the development of a mathematical modeling that can accurately replicate the interaction of water waves with porous coastal structures. Over the last 60 years, multiple techniques and solutions have been identiﬁed, from linearized solutions based on wave theories and constant friction coe ﬃ cients to very sophisticated Eulerian or Lagrangian solvers of the Navier-Stokes (NS) equations. In order to explore the ﬂow ﬁeld interior and exterior of the porous media under di ﬀ erent working conditions, the Smooth Particle Hydrodynamics (SPH) numerical simulation method was used to simulate the ﬂow distribution inside and outside a porous media applied to interact with the wave propagation. The ﬂow behavior is described avoiding Euler’s description of the interface problem between the Euler mesh and the material selected. Considering the velocity boundary conditions and the cyclical circulation boundary conditions at the junction of the porous media and the water ﬂow, the SPH numerical simulation is used to analyze the ﬂow ﬁeld characteristics, as well as the longitudinal and vertical velocity distribution of the back vortex ﬂow ﬁeld and the law of eddy current motion. This study provides innovative insights on the mathematical modelling of the interaction between porous structures and ﬂow propagation. Furthermore, there is a good agreement (within 10%) between the numerical results and the experimental ones collected for scenarios with porosity of 0.349 and 0.475, demonstrating that SPH can simulate the ﬂow patterns of the porous media, the ﬂow through the inner and outer areas of the porous media, and the ﬂow ﬁeld of the back vortex region. Results obtained and the new mathematical approach used can help to e ﬀ ectively simulate with high-precision the changes along the water depth, for a better design of marine and ocean engineering solutions adopted to protect coastal areas.


Introduction
In ocean and water conservancy projects, pore structures are widely used, such as reservoir dams, artificial reefs, and breakwaters. When water flows through pore structures, factors such as porosity rate and volume of porous media can be the major cause for different flow field fluctuations [1]. Therefore, to investigate the influence of porosity rate and volume of porous media on the flow field, which is crucial to characterize the features of flow field characteristics around the porous media [2], this paper analyzes the internal and external flow field characteristics of typical structures by simulating these physical phenomena.
The commonly used simulation method to study these effects is CFD (Computational Fluid Dynamics). CFD simulation can effectively simulate the flow field characteristics, however multiple

Smooth Particle Hydrodynamics (SPH) Methodology
The Smooth Particle Hydrodynamics method [12][13][14][15][16][17][18][19][20][21][22] is commonly used to describe a continuous fluid (or solid) as a group of interacting particles, despite each of these particles carries independent physical features such as different mass and dissimilar velocity. This method tends to characterize the continuous fluid by solving the behavior of the entire group of particles to simplify the problem and to determine the mechanical behavior of the entire system. In principle, as long as the number of particles is sufficient, the mechanical process can be accurately described. In practice, these particles are distributed irregularly, but in theory, they all move according to the motion law of the conservation governing equation. Any field of particles can be represented by applying a sufficiently accurate integral obtained by approximating the interpolation between adjacent particles. Considering this aspect, spatial derivatives such as gradient and scatter operators can be similarly estimated using the properties of all the particles in the Navier-Stokes equations, and the kernel approximation is commonly used to approximate the field function.
The procedure applied for this research uses the Smooth Particle Fluid Dynamics (SPH) method, which can be approximated by integral interpolation: where A(r) represent any function. For numerical simulations, the integral interpolation is calculated as follows: where a and b represent particles; m b and ρb represent particle mass and density, W ab is a kernel function, r is the distance between the reference particle and the adjacent (m), h is the smoothing length (m). To ensure a linear and angular momentum conservation, the asymmetric pressure gradient terms can be obtained as follows: where ∇ a W ab represents the kernel gradient taken relative to the position of the particle a and P represents the pressure; similarly, the divergence of the vector at a given particle a can be estimated by u: The viscous force can be calculated by combining the gradient and the divergence term, as follows: · u ab r ab · ∇ a W ab |r ab | 2 + η 2 (5) where υ represents the viscosity coefficient. Since pressure is highly sensitive to turbulence, it has been necessary to adopt the Laplace term in the pressure Poisson equation, written as follows [23]: P ab = P a − P b r ab = r a − r b The kernel function uses the cubic spline kernel function developed by Monaghan [24], which is related to the dimension of the field investigated and the distance between the relevant points. The kernel function is displayed below: where r is the distance between the reference particle and the adjacent one; q = r/h is the relative distance.

Equations for Flow Field Porous Media
Based on the previous studies conducted [19,25], the flow outside of the porous media is typically considered laminar and can be solved by the two-dimensional Navier-Stokes equations. The external flow field of pore structure is as follows: where ρ is the flow particle density; t is the time; u w is the external flow rate of the pore. Equation (10) is the governing equation of the external flow field of the porous media. This equation incorporates the influence of the current on the porous media during the movement, and the parameters to calculate the motion of the eddy current.
The internal flow field inside pore structure can be obtained as follows [19,25]: where u p1 is the conveying speed (m/s); n w is the porosity (/); K p is the permeability (/); C f is the nonlinear resistance coefficient (/).

Numerical Model Solving Process
In the first prediction step, the velocity and the position of the particles are calculated according to the momentum equations, and these equations can be expressed as follows [11]: r * = r t + u * ∆t (15) where ∆u is the particle velocity that changes during the prediction step; ∆t is the time increment; u t and r t represents the particle velocity and position at time t, respectively. The pressure term is based on the classical pressure Poisson equation that can be expressed as follows [23]: where ρ 0 represents the initial constant particle density; ρ * represents the central particle density after the prediction step, and P t+1 is the pressure of the particles at the t+1 time step. In the calibration of the second step, the pressure gradient term is combined with the momentum equation to ensure incompressibility. Pressure can be used to correct particle velocity as follows: where u t+1 is represents the particle velocity and r t+1 represents position at the moment of t + 1.

Free Surface Boundary
Shao and Lo [23] proposed a free surface treatment method where if the density of a particle is 10% lower than the reference density, it can be considered as a surface particle, and then zero pressure can be applied as a known boundary condition [8].
These assumptions are adopted for this study.

Fluid-Structure Coupling Boundary
The interface region is placed between the exchange flow and the porous flow region at the boundary of the pure fluid interface, and the interface line is located at the center of the region. Therefore, the width of the two regions (the fluid portion and the porous portion) at the interface boundary line is twice the distance between particles. By using the proposed SPH interface model, the average value of the interface region is obtained by solving the pressure Poisson equation in each calculation time step using the SPH kernel function. Only the adjacent particles located in the interface area are included in the summation.
At the same time, the interface boundary conditions can automatically satisfy the normal and tangential continuity.

Impermeable/Fixed Solid Wall Boundary
The impermeable/fixed solid wall is maintained along the boundary line. The pressure gradient between particles is adopted and when fixed virtual particles are used, the anti-skid boundary is selected.

Periodic Inflow and Outflow Boundaries Accompanied with a Damping Zone
According to previous studies conducted on pore structure logistics field, this paper adopts the boundary of inflow and outflow as displayed in Figure 1 below. between particles is adopted and when fixed virtual particles are used, the anti-skid boundary is selected.

Periodic Inflow and Outflow Boundaries Accompanied with a Damping Zone
According to previous studies conducted on pore structure logistics field, this paper adopts the boundary of inflow and outflow as displayed in Figure 1 below. The fluid particles near the left end of the channel can also interact with the fluid particles near the right end in the core affected area. In this sense, the upstream and downstream boundaries can actually be considered as overlapping. In the SPH numerical algorithm, we used the particle linked list to search for neighborhoods through a square grid with a side length of 2*H. Therefore, the fluid particles in the first column of the block located near the upstream boundary can be in contact with the particles in the last column of the block located near the downstream boundary, and vice versa, so the exit boundary flow velocity is also consistent with the inlet boundary flow velocity trend. As a preliminary test to verify this hypothesis, we set the calculation area as 3 m, the inlet flow rate as 0.00936 m 3 /s, and the particle size D as 0.1 m. In order to further test the performance of our method, three time stages t = 18.0 s, t = 20.0 s, and t = 22.0 s were considered to analyze its inlet and outlet flow rate, respectively. x = 0.0 represents the inlet flow rate distribution curve; x = 60.0 represents the outlet flow rate distribution curve. Analyzing the tests shown in Figure 2, it is possible to notice that the flow distribution of the particles is very similar between inlet (x = 0.0) and outlet (x = 60.0). Furthermore, the velocity ranges from 0 m/s to 1.5 m/s for each time step displayed proving the periodic boundary conditions of the inlet and outlet.  The fluid particles near the left end of the channel can also interact with the fluid particles near the right end in the core affected area. In this sense, the upstream and downstream boundaries can actually be considered as overlapping. In the SPH numerical algorithm, we used the particle linked list to search for neighborhoods through a square grid with a side length of 2*H. Therefore, the fluid particles in the first column of the block located near the upstream boundary can be in contact with the particles in the last column of the block located near the downstream boundary, and vice versa, so the exit boundary flow velocity is also consistent with the inlet boundary flow velocity trend. As a preliminary test to verify this hypothesis, we set the calculation area as 3 m, the inlet flow rate as 0.00936 m 3 /s, and the particle size D as 0.1 m. In order to further test the performance of our method, three time stages t = 18.0 s, t = 20.0 s, and t = 22.0 s were considered to analyze its inlet and outlet flow rate, respectively. x = 0.0 represents the inlet flow rate distribution curve; x = 60.0 represents the outlet flow rate distribution curve. Analyzing the tests shown in Figure 2, it is possible to notice that the flow distribution of the particles is very similar between inlet (x = 0.0) and outlet (x = 60.0). Furthermore, the velocity ranges from 0 m/s to 1.5 m/s for each time step displayed proving the periodic boundary conditions of the inlet and outlet. selected.

Periodic Inflow and Outflow Boundaries Accompanied with a Damping Zone
According to previous studies conducted on pore structure logistics field, this paper adopts the boundary of inflow and outflow as displayed in Figure 1 below. The fluid particles near the left end of the channel can also interact with the fluid particles near the right end in the core affected area. In this sense, the upstream and downstream boundaries can actually be considered as overlapping. In the SPH numerical algorithm, we used the particle linked list to search for neighborhoods through a square grid with a side length of 2*H. Therefore, the fluid particles in the first column of the block located near the upstream boundary can be in contact with the particles in the last column of the block located near the downstream boundary, and vice versa, so the exit boundary flow velocity is also consistent with the inlet boundary flow velocity trend. As a preliminary test to verify this hypothesis, we set the calculation area as 3 m, the inlet flow rate as 0.00936 m 3 /s, and the particle size D as 0.1 m. In order to further test the performance of our method, three time stages t = 18.0 s, t = 20.0 s, and t = 22.0 s were considered to analyze its inlet and outlet flow rate, respectively. x = 0.0 represents the inlet flow rate distribution curve; x = 60.0 represents the outlet flow rate distribution curve. Analyzing the tests shown in Figure 2, it is possible to notice that the flow distribution of the particles is very similar between inlet (x = 0.0) and outlet (x = 60.0). Furthermore, the velocity ranges from 0 m/s to 1.5 m/s for each time step displayed proving the periodic boundary conditions of the inlet and outlet.  It can also be observed from Figure 2 that the inlet and the exit boundary particles are basically consistent in the horizontal direction despite different time steps considered, confirming that the particles at the exit layer can smoothly enter the inlet across the boundary.
However, if the flow is significantly disturbed before exiting the outlet boundary, the inflow condition will inherit these disturbances and affect the downstream flow. In that case, at least a damping zone is needed downstream the inlet boundary and a numerical calibration is needed to Water 2020, 12, 918 7 of 25 clarify the effectiveness of the damping zone [26,27]. For this study, a damping zone is set at the junction of the exit boundary and the entrance boundary. Figure 3 displays the distribution of the velocity of the flow at the exit boundary along the depth of the water after passing through the damping zone. It can be seen that the speed variations remain relatively stable along the depth of the water at different times (t = 60 s; t = 70 s; t = 80 s; t = 90 s), hence the flow velocity at the outlet remains basically identical (the average flow velocity is stable at about 0.13 m/s at all times) as well as the inlet velocity. Under the effect of the periodic boundary, the inlet velocity can be stably maintained at 0.13 m/s so that the flow is significantly disturbed before it departures the exit boundary without affecting the downstream flow. This further validates the effectiveness of the damping region.
It can also be observed from Figure 2 that the inlet and the exit boundary particles are basically consistent in the horizontal direction despite different time steps considered, confirming that the particles at the exit layer can smoothly enter the inlet across the boundary.
However, if the flow is significantly disturbed before exiting the outlet boundary, the inflow condition will inherit these disturbances and affect the downstream flow. In that case, at least a damping zone is needed downstream the inlet boundary and a numerical calibration is needed to clarify the effectiveness of the damping zone [26,27]. For this study, a damping zone is set at the junction of the exit boundary and the entrance boundary. Figure 3 displays the distribution of the velocity of the flow at the exit boundary along the depth of the water after passing through the damping zone. It can be seen that the speed variations remain relatively stable along the depth of the water at different times (t = 60 s; t = 70 s; t = 80 s; t = 90 s), hence the flow velocity at the outlet remains basically identical (the average flow velocity is stable at about 0.13 m/s at all times) as well as the inlet velocity. Under the effect of the periodic boundary, the inlet velocity can be stably maintained at 0.13 m/s so that the flow is significantly disturbed before it departures the exit boundary without affecting the downstream flow. This further validates the effectiveness of the damping region. For all these reasons, it is possible to confirm the numerical characteristics of the periodic cyclic boundary of the flow.

Model Verification
In order for SPH to work in numerical sink simulations, characterizing the correct flow into the outflow boundary is key to achieve a continuous steady flow cycle. In this study, the flow was simulated using periodic boundary conditions and a damping zone. The position and velocity of the particles were firstly initialized, simulating the state of the current at the beginning of the movement. At the same time, during the interaction between the current and the porous media, the flow velocity was considered to change with the variation of the porosity. According to the experiments in the literature, porosity values of 0.349 and 0.475 [7] are the most commonly used and hence have been considered in this work to verify the model adopted.
The numerical water tank used for verification is shown in Figure 4 below [7]. The numerical flume has a pore structure placed 2.915 m away from the inlet boundary. The numerical model tank has a length of 6.0 m, an inlet flow rate of 0.13 m/s, a pore structure long 0.15 m and elevated 0.075 For all these reasons, it is possible to confirm the numerical characteristics of the periodic cyclic boundary of the flow.

Model Verification
In order for SPH to work in numerical sink simulations, characterizing the correct flow into the outflow boundary is key to achieve a continuous steady flow cycle. In this study, the flow was simulated using periodic boundary conditions and a damping zone. The position and velocity of the particles were firstly initialized, simulating the state of the current at the beginning of the movement. At the same time, during the interaction between the current and the porous media, the flow velocity was considered to change with the variation of the porosity. According to the experiments in the literature, porosity values of 0.349 and 0.475 [7] are the most commonly used and hence have been considered in this work to verify the model adopted.
The numerical water tank used for verification is shown in Figure 4 below [7]. The numerical flume has a pore structure placed 2.915 m away from the inlet boundary. The numerical model tank has a length of 6.0 m, an inlet flow rate of 0.13 m/s, a pore structure long 0.15 m and elevated 0.075 m, a sink channel slope of 0.005, and a viscosity of 1.0 × 10 −6 (Pa·s). Finally, the flow field characteristics were simulated at the porosity of 0.349 and 0.475.  It can be seen from Figure 5 that the flow simulation based on the periodic boundary basically simulates the correct natural movement of the water flow. The internal and external flow fields and back vortex flow fields of the porous media are also visible in Figure 4. The flow field characteristics, the flow velocity U, V, and water depth H were processed dimensionless to verify the flow field velocity distribution characteristics of the porous media at 0.475 and 0.349 and results are shown in Figures 6 and 7, respectively.    X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system.
The numerical simulation of the longitudinal and vertical velocity distributions at different locations around the porous media fit well (±5-10%) with the velocity curves measured for the experimental experiments. This confirms that the SPH method is feasible for the simulation of the flow field around the porous media and the analysis of the longitudinal and vertical velocity at a porosity of 0.475.
According to the above model verification displayed in Figure 7, it can also be seen that the SPH method can simulate the flow field inside and around the porous media at a porosity of 0.349, and results are consistent with the experiment of the flow field inside and outside the porous media. The model error is between 10% and 15%.  [7]. H: water depth (m); U: longitudinal flow rate (m/s), V: vertical flow rate (m/s), U0: Inlet boundary velocity (m/s). X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system. The numerical simulation of the longitudinal and vertical velocity distributions at different locations around the porous media fit well (±5-10%) with the velocity curves measured for the experimental experiments. This confirms that the SPH method is feasible for the simulation of the flow field around the porous media and the analysis of the longitudinal and vertical velocity at a porosity of 0.475.
According to the above model verification displayed in Figure 7, it can also be seen that the SPH method can simulate the flow field inside and around the porous media at a porosity of 0.349, and results are consistent with the experiment of the flow field inside and outside the porous media. The model error is between 10% and 15%.
Observing Figures 6 and 7, which show the average flow direction and vertical velocity components at different locations, it is also possible to notice specific features that are worth to be highlighted. At X/H = 0.33, the velocity gradient is large. In the vertical position, the upstream edge of the structure (X/H = 0.33) is 1 < Y/h < 1.5, and the velocity is small. At the downstream edge of the structure (X/H = 1.67), according to its velocity curve, it can be seen that the permeability of porous media is good, and its reverse velocity is small. X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system.
Observing Figures 6 and 7, which show the average flow direction and vertical velocity components at different locations, it is also possible to notice specific features that are worth to be highlighted. At X/H = 0.33, the velocity gradient is large. In the vertical position, the upstream edge of the structure (X/H = 0.33) is 1 < Y/h < 1.5, and the velocity is small. At the downstream edge of the structure (X/H = 1.67), according to its velocity curve, it can be seen that the permeability of porous media is good, and its reverse velocity is small.
X/H = 0.33 indicates that the position is 0.02475 m from the left boundary of the pore structure which is placed at the origin. This location is characterized by the velocity distribution of the flow field when the water flow is coupled with the left boundary of the pore structure. Because this is a periodic boundary and a damping region is added to the flow field, the damping region typically causes some disturbance to the flow field in the initial stage, but the effect of the damping region on the flow field velocity gradually decrease after the particles run for 20 s. It was verified within the model that the deviation has a small effect on the flow field around the pores, and it can be seen that the flow field around the pores is not greatly affected by the damping area.
However, in order to ensure that in the oncoming flow field of the pore structure, the flow simulation of the three parts of the pore flow field and the back vortex field of the pore structure have good convergence, the convergence analysis ( Figure 8) was completed. vertical flow rate (m/s), U0: Inlet boundary velocity (m/s). X is the distance between the particle and the left boundary of the pore structure, H is the height of the pore structure, and X/H is a dimensionless treatment of the X axis. The pore structure is placed at the origin of this system. X/H = 0.33 indicates that the position is 0.02475 m from the left boundary of the pore structure which is placed at the origin. This location is characterized by the velocity distribution of the flow field when the water flow is coupled with the left boundary of the pore structure. Because this is a periodic boundary and a damping region is added to the flow field, the damping region typically causes some disturbance to the flow field in the initial stage, but the effect of the damping region on the flow field velocity gradually decrease after the particles run for 20 s. It was verified within the model that the deviation has a small effect on the flow field around the pores, and it can be seen that the flow field around the pores is not greatly affected by the damping area.
However, in order to ensure that in the oncoming flow field of the pore structure, the flow simulation of the three parts of the pore flow field and the back vortex field of the pore structure have good convergence, the convergence analysis ( Figure 8) was completed.   Figure 8, it is possible to state that the particles pass through these three areas, and that due to the good convergence and good curve fitting that can be observed, the simulation with different particle sizes still meets the requirements of good convergence.
In practical engineering, the seepage problem of permeable structure will be helpful to the analysis of important engineering structures. It is necessary to study the flow state of fluid in permeable structure. The traditional velocity measurement method is difficult to be used in permeable buildings. Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are other robust experimental techniques to solve this problem [28][29][30][31][32][33][34][35]. In this study, the Lagrangian particle tracking method is used to measure the velocity distribution along the water depth at different positions.

Model Application
According to the model verification previously described in Section 4, it can be confirmed that the numerical simulation method based on SPH simulates appropriately the water flowing on the porous media, thereby it was decided to apply this model to a water tank and a porous media featuring periodic boundaries and a damping zone. The length of the domain is 60 m, at a porosity value is 0.5, particle size is 0.15, pore length, width, and height ranges within 0.6 m and 1.5 m, water  Figure 8 depicts the convergence analysis of the vortex area behind the pore structure and the vortex structure in the oncoming flow area when the porosity is 0.475 and 0.349, respectively. Analyzing Figure 8, it is possible to state that the particles pass through these three areas, and that due to the good convergence and good curve fitting that can be observed, the simulation with different particle sizes still meets the requirements of good convergence.
In practical engineering, the seepage problem of permeable structure will be helpful to the analysis of important engineering structures. It is necessary to study the flow state of fluid in permeable structure. The traditional velocity measurement method is difficult to be used in permeable buildings. Particle image velocimetry (PIV) and particle tracking velocimetry (PTV) are other robust experimental techniques to solve this problem [28][29][30][31][32][33][34][35]. In this study, the Lagrangian particle tracking method is used to measure the velocity distribution along the water depth at different positions.

Model Application
According to the model verification previously described in Section 4, it can be confirmed that the numerical simulation method based on SPH simulates appropriately the water flowing on the porous media, thereby it was decided to apply this model to a water tank and a porous media featuring periodic boundaries and a damping zone. The length of the domain is 60 m, at a porosity value is 0.5, particle size is 0.15, pore length, width, and height ranges within 0.6 m and 1.5 m, water depth considered corresponds to 2.4 m, viscosity is 1.0 × 10 −6 (Pa·s), structure coordinates are 30.0, 0.0. The schematic diagram of this configuration is displayed in Figure 9.
Water 2020, 12, x FOR PEER REVIEW 12 of 24 depth considered corresponds to 2.4 m, viscosity is 1.0 × 10 −6 (Pa·s), structure coordinates are 30.0, 0.0. The schematic diagram of this configuration is displayed in Figure 9.  Figure 10 shows the flow field velocity distribution of different volumes of porous media at the porosity of 0.5. According to the flow field displayed below, as the volume of the porous media increases (from top to bottom), the longitudinal flow velocity U of the upper layer of the porous media gradually increases, and the upward-incidence flow region is also expanding, but the enlarged region is not well defined. During the expansion process, the back vortex area expands, and multiple vortex motions are formed. The vortex area increases with the volume of the porous media, and begins to form when the vertical section of the porous media is equal to 1.2 m × 1.2 m. While there is an increase of pore structure volume, the back vortex area expands as well but the eddy strength decreases.  Figure 10 shows the flow field velocity distribution of different volumes of porous media at the porosity of 0.5. According to the flow field displayed below, as the volume of the porous media increases (from top to bottom), the longitudinal flow velocity U of the upper layer of the porous media gradually increases, and the upward-incidence flow region is also expanding, but the enlarged region is not well defined. During the expansion process, the back vortex area expands, and multiple vortex motions are formed. The vortex area increases with the volume of the porous media, and begins to form when the vertical section of the porous media is equal to 1.2 m × 1.2 m. While there is an increase of pore structure volume, the back vortex area expands as well but the eddy strength decreases.

Flow Field Velocity Distribution Diagram of Different Volumes of Porous Media
Water 2020, 12, x FOR PEER REVIEW 12 of 24 depth considered corresponds to 2.4 m, viscosity is 1.0 × 10 −6 (Pa·s), structure coordinates are 30.0, 0.0. The schematic diagram of this configuration is displayed in Figure 9.  Figure 10 shows the flow field velocity distribution of different volumes of porous media at the porosity of 0.5. According to the flow field displayed below, as the volume of the porous media increases (from top to bottom), the longitudinal flow velocity U of the upper layer of the porous media gradually increases, and the upward-incidence flow region is also expanding, but the enlarged region is not well defined. During the expansion process, the back vortex area expands, and multiple vortex motions are formed. The vortex area increases with the volume of the porous media, and begins to form when the vertical section of the porous media is equal to 1.2 m × 1.2 m. While there is an increase of pore structure volume, the back vortex area expands as well but the eddy strength decreases.  Figure 11 displays the horizontal and vertical velocity of the flow in the upstream area of 0.6 m × 0.6 m pore structure (the area where water flows through the pore structure and the vortex area). In this Figure 11, as well as for Figures 12-14 Figure 11 displays the horizontal and vertical velocity of the flow in the upstream area of 0.6 m × 0.6 m pore structure (the area where water flows through the pore structure and the vortex area). In this Figure 11, as well as for Figures 12-14   Moving from x = 0.0 m towards the pore structure (located at 30.0 m) when the distance between the water particles and the pore material becomes smaller, the longitudinal velocity u on the upper side of the pore material increases, as well as the internal velocity of the pore material. The average velocity of this area tends to be stable at about 1 m/s, and the vertical average velocity V ranges between −0.2-0.1 m/s. When x = 30.1-30.5 m, the water flows through the pore area. Increasing x, at the initial part or the pore structure, the longitudinal velocity u increases from 0.15 m/s to 1.5 m/s, showing the trend of short-term reflux. Therefore, it can be considered that the pore structure of 0.6 m × 0.6 m has an obvious energy dissipation effect, while in the upper water area of the pore structure, the flow velocity shows a slow upward trend and tends to be stable at 1.3 m/s. In the meantime, the vertical average velocity V ranges between −0.2 and 0.4 m/s. When x = 31.0-36.0 m, this stage is the eddy current area at the back of the pore structure. It can be seen from the change of the back vortex flow field that with the increase of the volume of porous media, there is a turning point at the speed of H (m) = 0.6 m. There are two stages of velocity curves, before and after the turning point, that can be described independently. The first stage represents the velocity curve inside the pore structure, and the second stage corresponds to the velocity curve outside the pore structure. When H is less than the height of pore structure, the velocity of the flow particles increases with the increase of distance from the right boundary of pore structure, with an average velocity of 0.5 m/s; when the water level H is greater than the height of the pore structure, the velocity decreases getting far away from the right boundary of the pore structure, while when x > 30.6 m, and H < 0.6, the velocity of the flow particles on the right boundary of pore structure decreases with the increase of distance from the right boundary of the pore structure. When x = 36.0 m, the flow particle velocity reaches the peak of 1.3 m/s. For this scenario, the vertical average velocity V is between −0.2 and 0.1 m/s. Figure 12 displays the horizontal and vertical velocity of the flow in the upstream area of 0.9 m × 0.9 m pore structure (the area where water flows through the pore structure and the vortex area).

Analysis of Flow Field Inside and Outside the Porous Media
Water 2020, 12, x FOR PEER REVIEW 14 of 24 media, there is a turning point at the speed of H (m) = 0.6 m. There are two stages of velocity curves, before and after the turning point, that can be described independently. The first stage represents the velocity curve inside the pore structure, and the second stage corresponds to the velocity curve outside the pore structure. When H is less than the height of pore structure, the velocity of the flow particles increases with the increase of distance from the right boundary of pore structure, with an average velocity of 0.5 m/s; when the water level H is greater than the height of the pore structure, the velocity decreases getting far away from the right boundary of the pore structure, while when x > 30.6 m, and H < 0.6, the velocity of the flow particles on the right boundary of pore structure decreases with the increase of distance from the right boundary of the pore structure. When x = 36.0 m, the flow particle velocity reaches the peak of 1.3 m/s. For this scenario, the vertical average velocity V is between −0.2 and 0.1m/s. Figure 12 displays the horizontal and vertical velocity of the flow in the upstream area of 0.9 m × 0.9 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25-29.0 m, with the increase of x, that is, when the distance between water particles and pore structure becomes smaller, the particle flow velocity in the upstream area of pore structure slightly decreases when the longitudinal velocity u on the upper side of pore structure is 0.6 m higher than the one recorded on the side, while the internal velocity growth of the pore structure is still almost negligible. The pore structure is more stable when the overall average velocity is 0.6 m higher than the side length, and the average up-flow velocity tends to be stable when it is close to 1 m/s, while the vertical average velocity V ranges between −0.2 and 0.1 m/s. With the increase of x, the longitudinal velocity u on the upper side of the porous material increases slightly, and the decrease observed on the left side is equal to the height of the porous material. The total velocity of the upwelling increases gradually, and tends to be stable near 0.9 m/s. The vertical velocity of the upwelling area remains almost stable at 0 m/s. When x = 25-29.0 m, with the increase of x, that is, when the distance between water particles and pore structure becomes smaller, the particle flow velocity in the upstream area of pore structure slightly decreases when the longitudinal velocity u on the upper side of pore structure is 0.6 m higher than the one recorded on the side, while the internal velocity growth of the pore structure is still almost negligible. The pore structure is more stable when the overall average velocity is 0.6 m higher than the side length, and the average up-flow velocity tends to be stable when it is close to 1 m/s, while the vertical average velocity V ranges between −0.2 and 0.1 m/s. With the increase of x, the longitudinal velocity u on the upper side of the porous material increases slightly, and the decrease observed on the left side is equal to the height of the porous material. The total velocity of the upwelling increases gradually, and tends to be stable near 0.9 m/s. The vertical velocity of the upwelling area remains almost stable at 0 m/s. When x = 30.1-30.5 m, water flows through the pore structure area. With the increase of x, the longitudinal velocity u of the fluid in the pore structure decreases gradually, but this phenomenon is not significant, and the reflux trend weakens. Therefore, it can be judged that the energy dissipation effect of the pore structure is slightly stronger when the side length is 0.9 m than when the side length is 0.6 m, while in the upper water area of the pore structure, the flow velocity shows a slow upward trend. When the height h is constant, the flow velocity decreases with the increase of x, and the overall flow velocity still increases with the increase of H. When x = 31.0-36.0 m, an eddy current area at the back of pore structure can be noticed. It can be seen from the change of the back vortex flow field that there is a turning point in the velocity distribution curve at H = 1.5 m. When H is less than the height of porous material, the velocity increases with the increase of x, and the average velocity is 0.25 m/s; when h is greater than the height of porous material, the velocity decreases with the increase of x, and reaches 1.7 m/s at x = 36.0 m, and the back vortex strength decreases with the increase of the distance from the hole. The vertical velocity is relatively stable, in the range of −0.15 and 0.1 m/s. Figure 13 displays the horizontal and vertical velocity of the flow in the upstream area of 1.2 m × 1.2 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25-29.0 m, with the increase of x, the particle velocity in the upstream area of the pore structure decreases slightly when the longitudinal velocity u on the upper side of the pore structure is 0.9 m longer than that on the side, and the average velocity of upwelling tends to be stable at the place close to 0.8 m/s, and the vertical average velocity V is between −0.15-0.1 m/s. When x = 30.1-30.5 m, water flows through the pore area. When x = 30.1-30.5 m, water flows through the pore structure area. With the increase of x, the longitudinal velocity u of the fluid in the pore structure decreases gradually, but this phenomenon is not significant, and the reflux trend weakens. Therefore, it can be judged that the energy dissipation effect of the pore structure is slightly stronger when the side length is 0.9 m than when the side length is 0.6 m, while in the upper water area of the pore structure, the flow velocity shows a slow upward trend. When the height h is constant, the flow velocity decreases with the increase of x, and the overall flow velocity still increases with the increase of H. When x = 31.0-36.0 m, an eddy current area at the back of pore structure can be noticed. It can be seen from the change of the back vortex flow field that there is a turning point in the velocity distribution curve at H = 1.5 m. When H is less than the height of porous material, the velocity increases with the increase of x, and the average velocity is 0.25 m/s; when h is greater than the height of porous material, the velocity decreases with the increase of x, and reaches 1.7 m/s at x = 36.0 m, and the back vortex strength decreases with the increase of the distance from the hole. The vertical velocity is relatively stable, in the range of −0.15 and 0.1 m/s. Figure 13 displays the horizontal and vertical velocity of the flow in the upstream area of 1.2 m × 1.2 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25-29.0 m, with the increase of x, the particle velocity in the upstream area of the pore structure decreases slightly when the longitudinal velocity u on the upper side of the pore structure is 0.9 m longer than that on the side, and the average velocity of upwelling tends to be stable at the place close to 0.8 m/s, and the vertical average velocity V is between −0.15-0.1 m/s. When x = 30.1-30.5 m, water flows through the pore area. With the increase of x, the longitudinal velocity U of the fluid in the pore structure is gradually stable at 0.1 m/s, and the reflux trend is not strong. It can be seen that for 1.2 m × 1.2 m porous material, With the increase of x, the longitudinal velocity U of the fluid in the pore structure is gradually stable at 0.1 m/s, and the reflux trend is not strong. It can be seen that for 1.2 m × 1.2 m porous material, the energy dissipation effect of porous material is also stronger than before. When x = 31.0-36.0 m, when H is greater than the hole height, the velocity decreases with the increase of x, and the velocity at the right boundary of the pore structure decreases with the increase of the distance from the right boundary of the hole to 2.0 m/s. The back vortex intensity also decreases with the increase of x. The vertical velocity is relatively stable, in the range of −0.20 and 0.35 m/s. Finally, Figure 14   In summary, based on the changes of up-flow field tested, it can be stated that by increasing the pore volume, the up-flow velocity typically decreases and the energy dissipation becomes higher, but changes associated with the vertical velocity are not significant with the change of volume for the In summary, based on the changes of up-flow field tested, it can be stated that by increasing the pore volume, the up-flow velocity typically decreases and the energy dissipation becomes higher, but changes associated with the vertical velocity are not significant with the change of volume for the pore structure. Focusing on the change of the flow field inside and outside the hole it is possible to confirm that the internal velocity of the pore is mainly in the range of 0 and 0.2 m/s. With the increase of the porous media' volume, the internal velocity difference of the pore decreases gradually, while the particle velocity directly above the pore structure increases with the increase of the pore volume, and the velocity difference of the upper layer flow also gradually increases.
Furthermore, it can be seen from the change of the back vortex flow field that, with the increase of the volume of porous media, there are inflexion points at H (m) = 0.6, 0.9, 1.2 and 1.5 for the four groups of horizontal velocity, respectively, thus forming multiple velocity curves before and after the inflexion point. When H is less than the height of the hole, the flow velocity increases with the increase of x, however when H is greater than the height of the hole, the flow velocity increases with the increase of x. The intensity of the back vortex also decreases with the increase of x. From the change of vertical velocity, it can be observed that the velocity basically varies from −0.5 m/s to 0.5 m/s.

Longitudinal and Vertical Flow Field Distribution under Different Porosity
Porosity is also an important factor influencing the characteristics of the flow field in the porous media [36]. Therefore, the pore size was 1.2 × 1.2 (m), and the oncoming flow and reflux vortices were analyzed for porosity values of 0.5, 0.25 and 0.125, respectively. The velocity distribution and the flow field along the water depth were analyzed for each porosity investigated. Figure 15 displays the flow field in the upstream area of pore structure under different porosities 0.5, 0.25, and 0.125 for x = 25.0-29.0. As the porosity decreases, it can be seen that the longitudinal velocity curve of the upstream flow tends to be gradually stable.
Water 2020, 12, x FOR PEER REVIEW 17 of 24 pore structure. Focusing on the change of the flow field inside and outside the hole it is possible to confirm that the internal velocity of the pore is mainly in the range of 0 and 0.2 m/s. With the increase of the porous media' volume, the internal velocity difference of the pore decreases gradually, while the particle velocity directly above the pore structure increases with the increase of the pore volume, and the velocity difference of the upper layer flow also gradually increases. Furthermore, it can be seen from the change of the back vortex flow field that, with the increase of the volume of porous media, there are inflexion points at H (m) = 0.6, 0.9, 1.2 and 1.5 for the four groups of horizontal velocity, respectively, thus forming multiple velocity curves before and after the inflexion point. When H is less than the height of the hole, the flow velocity increases with the increase of x, however when H is greater than the height of the hole, the flow velocity increases with the increase of x. The intensity of the back vortex also decreases with the increase of x. From the change of vertical velocity, it can be observed that the velocity basically varies from −0.5 m/s to 0.5 m/s.

Longitudinal and Vertical Flow Field Distribution under Different Porosity
Porosity is also an important factor influencing the characteristics of the flow field in the porous media [36]. Therefore, the pore size was 1.2 × 1.2 (m), and the oncoming flow and reflux vortices were analyzed for porosity values of 0.5, 0.25 and 0.125, respectively. The velocity distribution and the flow field along the water depth were analyzed for each porosity investigated. Figure 15 displays the flow field in the upstream area of pore structure under different porosities 0.5, 0.25, and 0.125 for x = 25.0-29.0. As the porosity decreases, it can be seen that the longitudinal velocity curve of the upstream flow tends to be gradually stable.                    In summary, it can be seen from Figures 15-20 that as the porosity increases, the longitudinal flow velocity of the oncoming flow slightly changes, but the overall change is not significant, plus the longitudinal flow velocity of the back vortex decreases. Considering the vertical flow velocity of the upstream area, despite the increase of porosity, this parameter remains within the range between −0.2 m/s and 0.2 m/s, and the vertical flow velocity through the porous media decreases with the increase of porosity. From the vertical flow velocity of the back vortex, as the porosity increases, the vertical flow velocity tends to decrease slightly overall.
Due to the infiltration of the water, the permeable structure has a significant impact on its discharge capacity [37]. It can be seen from the horizontal and vertical comparison of the flow field in Figure 16 that with the increase of porosity, the flow resistance decreases. The distribution of longitudinal velocity of porous material changes in the upper layer of porous material.
The maximum turbulence intensity near the surface and behind the structure increases proportional to the porosity [38] at the same downstream position.

Convergence Verification of Pore Logistics Field Model
In order to verify the convergence of the SPH in the numerical model, the flow field characteristics of each section previously investigated were analyzed by applying multiple particle size distributions [8]. (a) (b) (c) Figure 22. Vertical velocity distribution of the inflow field with different particle spacing (a), vertical velocity distribution of the sea current flowing through the porous media with different particle spacing (b), and vertical velocity distribution of the rear vortex field with different particle spacing (c).
It was previously noticed that when the height of the pores is 1.2 m, the flow direction is opposite to the others when D = 0.096. This elevation corresponds to the upper flow area of the pore structure and justify the presence of turbulence as shown in the graph. However, results displayed in Figures   Figure 21. Longitudinal velocity distribution of the inflow field with different particle spacing (a), longitudinal velocity distribution of the sea current flowing through the porous media with different particle spacing (b), and longitudinal velocity distribution of the rear vortex field with different particle spacing (c). (a) (b) (c) Figure 22. Vertical velocity distribution of the inflow field with different particle spacing (a), vertical velocity distribution of the sea current flowing through the porous media with different particle spacing (b), and vertical velocity distribution of the rear vortex field with different particle spacing (c).
It was previously noticed that when the height of the pores is 1.2 m, the flow direction is opposite to the others when D = 0.096. This elevation corresponds to the upper flow area of the pore structure and justify the presence of turbulence as shown in the graph. However, results displayed in Figures   Figure 22. Vertical velocity distribution of the inflow field with different particle spacing (a), vertical velocity distribution of the sea current flowing through the porous media with different particle spacing (b), and vertical velocity distribution of the rear vortex field with different particle spacing (c).
It was previously noticed that when the height of the pores is 1.2 m, the flow direction is opposite to the others when D = 0.096. This elevation corresponds to the upper flow area of the pore structure and justify the presence of turbulence as shown in the graph. However, results displayed in Figures 21 and 22 show that the velocity distributions of the internal and external flow fields formed by the oncoming flow field, the back vortex flow field and the water flow through the porous media are consistent despite different particle spacing and the error range is between 10% and 15%, hence providing an acceptable grade of convergence.
Finally, Figures 23 and 24 show the back vortex field of the pore structure when the porosity is 0.349 and 0.475, respectively.  It can be seen from the flow field distribution diagrams in Figures 23 and 24 that the energy in the vicinity of the pore structure presents a dissipative trend and shows different vorticities under different porosities. When the porosity is 0.475, the energy dissipation in the vicinity of the pore structure is weak, the vortex flow pattern is more noticeable and multiple vortex flow patterns are formed on the back of the pore structure. However, when the porosity is 0.349, the energy dissipation near the pore structure is relatively strong, and the intensity of the vortex flow is relatively weak.
To enable academic colleagues to replicate the presented numerical data and have a better estimation of the numerical capabilities of the numerical method, details about the total number of  It can be seen from the flow field distribution diagrams in Figures 23 and 24 that the energy in the vicinity of the pore structure presents a dissipative trend and shows different vorticities under different porosities. When the porosity is 0.475, the energy dissipation in the vicinity of the pore structure is weak, the vortex flow pattern is more noticeable and multiple vortex flow patterns are formed on the back of the pore structure. However, when the porosity is 0.349, the energy dissipation near the pore structure is relatively strong, and the intensity of the vortex flow is relatively weak.
To enable academic colleagues to replicate the presented numerical data and have a better estimation of the numerical capabilities of the numerical method, details about the total number of It can be seen from the flow field distribution diagrams in Figures 23 and 24 that the energy in the vicinity of the pore structure presents a dissipative trend and shows different vorticities under different porosities. When the porosity is 0.475, the energy dissipation in the vicinity of the pore structure is weak, the vortex flow pattern is more noticeable and multiple vortex flow patterns are formed on the back of the pore structure. However, when the porosity is 0.349, the energy dissipation near the pore structure is relatively strong, and the intensity of the vortex flow is relatively weak.
To enable academic colleagues to replicate the presented numerical data and have a better estimation of the numerical capabilities of the numerical method, details about the total number of particles used, the particle size, the physical simulation times plus the CPU time and the CPU cores used are displayed and summarized in Table 1

Conclusions
Nowadays, a variety of man-made structures such as submerged breakwaters, fishing reefs, outfall protections, armor layers, rubble-mound or berm breakwaters, built of gravel or artificial units, can be found in coastal areas worldwide for multiple purposes (e.g., coastal flooding protection, coastal erosion protection). All these structures are characterized by the fact that either some of their layers or their full structure are porous, hence they are permeable to flows induced by waves and currents. Therefore, the hydrodynamics, the stability and performance of these structures are dependent on the characteristics of the waves and their interaction with permeable material and it is crucial to provide mathematical formulations to model these features and the role played by the porous material considering its geometry, its location, and its characteristics (which could influence the wave propagation, diffusion, overtopping, and the consequent turbulence generated). This paper enabled the calibration and validation of a SPH model utilized to characterize the influence of porous media on flow propagation and outcomes can be summarized as follows: 1. Based on the analysis of the area upstream the porous media, the longitudinal flow velocity of the ascending flow decreases slightly with the increase of the volume of the pore structure; 2. Based on the analysis of the internal and external flow fields of the porous media, the internal flow velocity of the porous media for the configurations tested is in the range of 0-0.2 m/s. It was possible to notice that with the increase of the volume of the porous media, the flow velocity of the upper layer of the porous media firstly tends to accelerate and then stabilizes. The longitudinal flow velocity of the upper layer of the flow increases with the increase of water depth, while the vertical velocity fluctuates sharply at the fluid-solid boundary; 3. Based on the back vortex analysis, there is an inflection point in the flow velocity distribution of the upper and lower layers of the porous media. As the distance from the porous media increases, the flow velocity increases, so that the velocity distribution curve of the back vortex flow forms different concavities and convexities before and after the porous media, and as the pore volume gradually increases, a plurality of vortices are formed on the backflow surface. With the rotary motion, the area of the back vortex increases as the volume of the pore increases, but the strength decreases as the volume increases. 4. Overall, with the increase of porosity, the longitudinal flow and back vortex flow slightly decrease in the vertical direction, and the SPH-based model calculation has good convergence (error between experimental and numerical results within 10%).
Funding: This research was funded by Zhejiang Province Natural Science Foundation, grant number LQ16E090001, the National Natural Science Foundation of China, grant number 51509134.