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

Abstract

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 identified, from linearized solutions based on wave theories and constant friction coefficients to very sophisticated Eulerian or Lagrangian solvers of the Navier-Stokes (NS) equations. In order to explore the flow field interior and exterior of the porous media under different working conditions, the Smooth Particle Hydrodynamics (SPH) numerical simulation method was used to simulate the flow distribution inside and outside a porous media applied to interact with the wave propagation. The flow 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 flow, the SPH numerical simulation is used to analyze the flow field characteristics, as well as the longitudinal and vertical velocity distribution of the back vortex flow field and the law of eddy current motion. This study provides innovative insights on the mathematical modelling of the interaction between porous structures and flow 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 flow patterns of the porous media, the flow through the inner and outer areas of the porous media, and the flow field of the back vortex region. Results obtained and the new mathematical approach used can help to effectively 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.

1. 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 complications and lengthy simulations are associated with the establishment of a proper mesh [3,4]. Boundary element simulations can effectively simulate the influence of water flow on pore structure, but these methods only divide the elements on the boundary domain to meet the approximations applied to the functions implemented for the boundary conditions of the control equation. Therefore, compared with the finite element method, the boundary element method has the advantage of fewer elements to consider and a simpler data preparation. However, when the boundary element method is used to solve the nonlinear problem, which corresponds to the integral corresponding to the nonlinear term, more complications are frequently encountered because this kind of integration has strong singularity near single points, making an accurate solution very challenging to be obtained. Therefore, to reduce and avoid these complications, the key properties of the Smooth Particle Hydrodynamics (SPH) method, which is based on the pure Lagrangian numerical simulation method, were adopted for this study. This method can track the trajectory of each particle, and has good advantages in terms of visualization, high precision, and large processing. Many studies have been completed to date to explore the hydrodynamic characteristics of porous media. For example, Hossei Basser et al. [5] used SPH method to characterize the gravity flow of a single fluid in porous media and multi-fluid flows in impermeable layers and porous media, as well as the influence of unidirectional flow on pore structure and particle motion. The results that were obtained confirmed how the SPH numerical simulation method had the characteristics of high accuracy and visualization, which could effectively and accurately simulate the movement of each particle in the flow field inside and outside the pore structure, describing very accurately the characteristics of the flow field inside and outside the pore structure.

Pore flow field simulation can be simulated focusing on two aspects: internal flow field and external flow field, which refer to the specific flow field formed in large-scale ocean engineering, such as artificial reefs, dam break and breakwater, under the condition of laminar and turbulent flows. Boundary conditions are very important and are typically periodic or open. The first category means that the flow particles run from the inlet boundary to the outlet boundary, re-entering the inlet boundary and form a continuous loop. The second category refers to the flow particles moving in from the inlet boundary to the outlet boundary without re-entering from the inlet boundary, because this approximation considers that new flow particles will move in following an uninterrupted process. Under the open boundary conditions, researchers typically simulate the internal and external flow fields, effectively realizing the continuous process of particles motion [6]. However, by using these conditions although the accuracy of the inlet and outlet flows is high, the interaction with the porous media is relatively complex and affected by lower accuracy. Therefore, periodic boundary conditions accompanied with a damping zone were employed in this paper to complete the work mentioned above.

When simulating pore structure models, porosity must be considered, because factors such as pore permeability can affect the flow field around the pore, like for example in processes such as dam break, where the dam body is composed of pore structure. Therefore, under different porosity, the dam break strength is different, generating different impacts in case of failure of the dam taken into consideration. Additionally, when the flow comes in contact with the pore structure, the seepage characteristics of the pore structure can further affect the surrounding flow field [7], so the porosity that affects the pore structure is one of the key factors to be analyzed.

To date, multiple studies have been conducted on the characterization of pore structures. Qinqin Gui et al. [8] studied the coupling law of wave and pore structure by using the incompressible SPH method. Through numerical simulation and experiment, the typical flow around and inside the underwater porous structure was compared, however the analysis of the internal and external flow field of the porous structure was not considered as well as the influence of water flow on the porous structure. Abbas Khayyer et al. [9] proposed a numerical solution for variable porosity, because the purpose of this study was to investigate the effect of waves on structures with different porosity. Vanneste and Troch [10] conducted a detailed two-dimension numerical study on the wave propagation and verified the calculation results in terms of average water level and spatial pore pressure distribution from the incident wave field to the normal wave field. Shao [11] proposed an Incompressible Smooth Particle Hydrodynamics (ISPH) method to simulate the interaction between waves and porous media, and achieved interesting results on the analysis of flows and pore structures: through comparative analysis, the error obtained by Shao [11] was between 10% and 15%, which better reflects the accuracy of numerical simulations when the SPH is applied [8]. To summarize what was found in literature, to date, traditional SPH method is generally used to estimate the influence of waves on the porous media, but seldom for the coupling the effects of water flow and pore separately. Therefore, there is a need for research to be conducted to demonstrate that the SPH method can be adopted to simulate the action process of water flow on porous materials and the consequent effects. To fulfill this gap, SPH is used in this work to simulate coupling effects of water flow and porous media with the periodic circulation flow as a boundary condition. The influence of the flow field on the porous media was analyzed, together with the flow field inside and outside the porous media, the longitudinal and vertical velocity distribution, and the eddy current law of the eddy current field.

2. 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:

$$A(r)={\displaystyle \int A(r^{\prime })}W(r-r^{\prime },h)d{r}^{\prime }$$

(1)

where A(r) represent any function. For numerical simulations, the integral interpolation is calculated as follows:

$$A_a(r)={\displaystyle \sum {}_b{m}_b}\frac{A_b(r)}{{\rho }_b}{W}_{ab}$$

(2)

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:

$${(\frac{1}{\rho }\nabla P)}_a={\displaystyle \sum {}_b{m}_b}(\frac{P_a}{P_{a^2}}+\frac{P_b}{P_{b^2}}){\nabla }_a{W}_{ab}$$

(3)

where ∇_aW_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:

$$\nabla u_a={\rho }_a{\displaystyle \sum {}_b{m}_b}(\frac{u_a}{P_{a^2}}+\frac{u_b}{P_{b^2}}){\nabla }_a{W}_{ab}$$

(4)

The viscous force can be calculated by combining the gradient and the divergence term, as follows:

$${({\upsilon }_0{\nabla }^{2}u)}_a={\displaystyle \sum {}_b{m}_b}\frac{2({\upsilon }_a+{\upsilon }_b)}{({\rho }_a+{\rho }_b)}\cdot \frac{u_{ab}{r}_{ab}\cdot {\nabla }_a{W}_{ab}}{{\left|r_{ab}\right|}^2+{\eta }^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]:

$$\nabla \cdot {(\frac{\nabla P}{\rho })}_a={\displaystyle \sum {}_b{m}_b}\frac{8\cdot P_{ab}{r}_{ab}\cdot {\nabla }_a{W}_{ab}}{{({\rho }_a+{\rho }_b)}^2\cdot ({\left|r_{ab}\right|}^2+{\eta }^2)}$$

(6)

$$P_{ab}=P_a-P_b\hspace{1em}\hspace{1em}{r}_{ab}=r_a-r_b$$

(7)

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:

$$W(r,h)=\{\begin{array}{ccc}(10\times (1-1.5q^2+0.75q^3))/7\pi h^2& & 0\le q< 1\\ (10{(2-q)}^3)/28\pi h^2& & 1\le q\le 2\\ 0& & q>2\end{array}$$

(8)

where r is the distance between the reference particle and the adjacent one; q = r/h is the relative distance.

3. SPH Model

3.1. 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:

$$\frac{d\rho }{\rho dt}+\nabla \cdot u_{\mathrm{w}}=0$$

(9)

$$\frac{d{u}_w}{dt}= -\frac{1}{\rho }\nabla \mathrm{P}+\mathrm{g}+ \mathsf{\upsilon }{\nabla }^2{u}_w$$

(10)

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]:

$$\frac{d{u}_p}{dt}=-\frac{1}{\rho }\nabla P+g+\upsilon \cdot {\nabla }^2{u}_p-\frac{C_f{n}_w{}^2}{\sqrt{K_P}}{u}_p\left|u_p\right|$$

(11)

$$K_p=1.643\times {10}^{-7}{\left[\frac{d_{50}}{d_0}\right]}^{1.57}\left[\frac{n_w{}^3}{{(1-n_w)}^2}\right] C_f=100{[d_{50}(m)\cdot {\left(\frac{n_w}{K_p}\right)}^{1/2}]}^{-1.5}$$

(12)

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 (/).

3.2. 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]:

$$\mathsf{\Delta }{u}_{*}={(g+\upsilon {\nabla }^2{u}_p-\frac{\upsilon n_w}{K_P}u-\frac{C_f{n}_w{}^2}{\sqrt{K_P}}u\left|u\right|)}_t\mathsf{\Delta }t$$

(13)

$$u_{*}=u_t+\mathsf{\Delta }{u}_{*}$$

(14)

$$r_{*}=r_t+u_{*}\mathsf{\Delta }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]:

$$\nabla \cdot (\frac{1}{\rho }\nabla P_{t+1})=\frac{{\rho }_0-{\rho }_{*}}{{\rho }_0\mathsf{\Delta }{t}^2}$$

(16)

where ρ₀ 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:

$$\mathsf{\Delta }{u}_{**}=\frac{-1}{{\rho }_{*}}\nabla p_{t+1}\mathsf{\Delta }t$$

(17)

$$u_{t+1}=u_{*}+\mathsf{\Delta }{u}_{**}$$

(18)

$$r_{t+1}=r_t+\frac{u_{t+1}+u_t}{2}\mathsf{\Delta }t$$

(19)

where u_t+1 is represents the particle velocity and r_t+1 represents position at the moment of t + 1.

3.3. Boundary Conditions

3.3.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.

3.3.2. 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.

3.3.3. 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.

3.3.4. 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.

Particles leaving one side of the computational domain move directly through the other side, and particles near the open side boundary interact with particles near the complementary open side boundary. 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³/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 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.

4. 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⁻⁶ (Pa·s). Finally, the flow field characteristics were simulated at the porosity of 0.349 and 0.475.

The flow fields of velocity vectors for porosity equal to 0.349 (Figure 5a) and 0.475 (Figure 5b) are shown below in Figure 5.

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 Figure 6 and Figure 7, respectively.

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 Figure 6 and Figure 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.

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.

5. 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⁻⁶ (Pa·s), structure coordinates are 30.0, 0.0. The schematic diagram of this configuration is displayed in Figure 9.

5.1. Flow Field Velocity Distribution Diagram of Different Volumes of Porous Media

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.

5.2. Analysis of Flow Field Inside and Outside the Porous Media

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 Figure 12, Figure 13 and Figure 14, x = 25–29.0 m represents the horizontal coordinate range of the upstream area, x = 30.1–30.5 m represents the horizontal coordinate range of the flow through the pore structure, and x = 31.0–36.0 m represents the horizontal coordinate range of the vortex area.

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).

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.

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 displays the horizontal and vertical velocity of the flow in the upstream area of 1.5 m × 1.5 m pore structure (the area where water flows through the pore structure and the vortex area). When x = 25–29.0 m, the overall velocity of the upwelling increases gradually, and tends to be stable (approximately 0.75 m/s), while the vertical velocity ranges between −0.2 and 0.2 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 tends to 0, and the reflux trend is not noticeable. It can be concluded that the hole has obvious energy dissipation effect on the hole of 1.5 m × 1.5 m, while in the upper water area of the hole, the velocity keeps rising slowly. With the increase of distance, the water content of the upper part of the pore structure shows a slowly rising trend, the vertical velocity increases at the same time. When x = 31.0–36.0 m, when h is less than the height of porous material, the flow velocity increases with the increase of the right boundary distance of porous material, and the average flow velocity is approximately 0.02 m/s; when h is greater than the height of porous material, the flow velocity decreases with the increase of the right boundary distance of porous material, and the flow velocity at the right boundary of porous material reaches 2.2 m/s, and the back vortex strength decreases with the increase of the distance from the hole. The vertical velocity is relatively stable, in the range between −0.35 and 0.35 m/s.

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.

5.3. 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.

Figure 16 shows the flow field when the water flows through the pore structure under different porosities 0.5, 0.25 and 0.125 for x = 30.1–30.5. The longitudinal velocity difference of the internal flow field of the pore structure increases with the decrease of porosity.

Figure 17 presents the flow field in the back vortex region with different porosities 0.5, 0.25 and 0.125 for x = 31.0–36.0. With the decrease of porosity, the flow rate changes are not significant, but the back vortex region increases slightly, and the flow rate at the inflection point also initially increases but later decreases in combination with the increase of porosity. When the porosity is 0.25, the particle longitudinal velocity at the inflection point is the maximum calculated.

Figure 18 displays the vertical velocity distribution of the flow field in the upstream area of the pore structure under different porosities 0.5, 0.25, and 0.125 for x = 25.0–29.0. With the gradual decrease of porosity, it can be seen that the vertical velocity difference of each position of the upstream flow is gradually reduced.

Figure 19 shows the vertical velocity distribution of the flow field inside and directly above the pore structure under different porosities 0.5, 0.25, and 0.125 for x = 30.1–30.5. With the gradual decrease of porosity, it can be seen that in the fluid solid interface area, the change of horizontal velocity becomes more evident, and the velocity reduces.

Figure 20 presents the vertical velocity distribution of the back vortex flow field under different porosities 0.5, 0.25 and 0.125 for x = 31.0–36.0. As the porosity decreases gradually, it can be seen that when the height H is constant, the vertical velocity difference at each position is characterized by a gradually increasing trend.

In summary, it can be seen from Figure 15, Figure 16, Figure 17, Figure 18, Figure 19 and Figure 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.

5.4. 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]. Up-flow area, pore flow area, back vortex area, and the horizontal and vertical velocity distribution along the depth were calculated with the particle size D of 0.096 m, 0.01 m, and 0.104 m as shown in Figure 21 and Figure 22 to verify convergence, at the selected positions x = 28.0 m, 30.3 m, 33 m.

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 Figure 21 and Figure 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, Figure 23 and Figure 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 Figure 23 and Figure 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. Simulation details for each configuration tested.

Particle Size D (m)	Total Number of Particles	Total Physical Simulation Time (s)	CPU Time	CPU Cores
0.0104	13,914	100	27 h 05 min	4
0.01	15,000	100	27 h 14 min	4
0.096	16,225	100	8 h 28 min	20
0.1	15,000	100	8 h 23 min	20
0.104	13,914	100	8 h 18 min	20

for all the configurations tested (D = 0.0104 m, D = 0.01 m, D = 0.096 m, D = 0.1 m, and D = 0.104 m).

6. 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%).