Fully Coupled Simulation of Interactions Among Waves, Permeable Breakwaters and Seabeds Based on N−S Equations

Interstitial flows in breakwater cores and seabeds are a key consideration in coastal and marine engineering designs and have a direct impact on their structural safety. In this paper, a unified fully coupled model for wave−permeable breakwater−porous seabed interactions is built based on an improved N−S equation. A numerical wave flume is constructed, and numerical studies are carried out by applying the finite difference method. In combination with a physical model test, the accuracy of the numerical simulation results is verified by comparing the calculated and measured values of wave height at measurement points and the seepage pressure within the breakwater and seabed. On this basis, the characteristics of the surrounding wave field and the internal flow field of the pore structure, as well as the evolution process of the fluctuating pore water pressure inside the breakwater and seabed, are further analyzed. The spatial distribution of the maximum fluctuating pore water pressure in the breakwater is compared between two cases by considering whether the seabed is permeable, and then the effect of seabed permeability on the dynamic pore water pressure in the breakwater is clarified. This study attempts to provide a reference for breakwater design and the protection of nearby seabeds.


Introduction
Permeable breakwaters, also known as mound rubble breakwaters, are important and commonly used coastal structures. In terms of the causes of structural instability, there are two forms of structural failure, namely, the partial instability of the structure itself under wave actions and the structural instability caused by wave-induced responses of the seabed in the vicinity of the structure. The stability of a permeable breakwater is closely related to the characteristics of the flow field inside and around it. When water waves act on these permeable structures, part of the energy enters the breakwater and creates seepage. If the seabed is also permeable, the presence of the seabed will affect the hydrodynamic properties within the breakwater, such as the distribution of flow velocity and pressure. The presence of a coastal structure, on the other hand, also affects both the surrounding wave field and the seabed seepage. The wave−permeable breakwater−seabed interaction is a triple coupled process, and its study is of great engineering importance.
At present, there are very limited experimental studies on the physical modeling of wave−seabed−structure interactions. Instead, numerical simulations are more frequently used, among which many scholars have adopted a partitioned modeling approach where pore flow and external flow fields are modeled separately, and then the data are transferred or coupled through boundary conditions, such as the combined N−S equation model and Biot's consolidation theory (Liu and Garcia, 2007) and the combined Reynolds Averaged Navier-Stokes model considering the turbulence effect and Biot's theory (Cheng et al., 2007;Ye et al., 2015;Zhao and Jeng, 2015). Different methods are usually adopted for different regions, and the finite element method is often used for assessing pore structure, such as the BEM-FEM (Mostafa et al., 1999) and the FVM-FEM (Zhao et al., 2018).
In the partition modeling strategy, it is difficult to realize a fully coupled simulation due to the unknown boundary conditions on the interface between the external wave motion and the pore water flow, so an iterative algorithm is required. Thus, weakly coupled models have been adopted in many studies because the response of the pore structure can be studied. Another type of approach models the interaction by adding flow resistance terms provided by the pore structure in a certain portion of the computational region. The governing equation can usually be expressed in a uniform form throughout the whole flow field. Sakakiyama and Kajima (1992) simulated wave−permeable breakwater interactions by adding nonlinear drag force terms of inertial and velocity squared resistance in the N−S equation and determined the resistance coefficient through the reflection and transmission coefficients of the breakwater to waves. Mizutani et al. (1996) also developed a pore flow model by the same method in which the combined effects of inertial force, laminar flow resistance and turbulent flow resistance were considered to study the interactions between waves and the breakwater, and the potential flow theory was applied for the external flow field, which was solved numerically by a combination of the BEM and FEM. It was found that the force of nonbreaking waves on the submerged breakwater armors was greater than that of the breaking waves. Wang (2012) modeled the interaction of a wave field and a permeable breakwater based on an improved N−S equation, solved it in a numerical wave flume (Kawasaki, 1999) using the finite difference method, and verified the computational results by comparison with physical model tests. Ren et al. (2016) also developed a Volume Averaged and Favre Averaged Navier-Stokes model of wave-permeable breakwater interaction by adding laminar flow resistance and turbulent flow resistance, applied an improved SPH method to solve it, and obtained good simulation results.
In the above research, only the interactions between waves and permeable breakwaters were considered. However, the permeability of the seabed sometimes is not negligible. Therefore, this paper further develops the N−S equation and establishes a unified governing equation on different flow regions to simulate a fully coupled wave−permeable breakwater−porous seabed interaction. To consider the impact of the permeability of breakwaters and seabeds, corresponding flow resistance terms are added to their respective regions, which are distinguished by a porosity parameter, and thus an explicit boundary condition between the wave motion and the porous flow is avoided. An explicit solution format of finite difference is applied to calculate the model. A corresponding physical model test is implemented to verify the calculation results. Finally, a detailed analysis of the flow field in and around the breakwater and seabed is presented, and the dynamic pressure distribution is discussed.

Governing equations
Suppose that the soil skeleton of the seabed and the breakwater are rigid porous media and that the fluid is viscous and incompressible. Using the Morison equation for reference, the effect of pore material on the flow is reflected in the motion equation in the form of a pressure difference that consists of two parts: the velocity term proportional to the square of fluid velocity and the inertia term proportional to the acceleration. In addition, Mizutani et al. (1996) also considered a linear drag force term similar to Darcy flow based on the Morison equation. The unified governing equations of the motion of external fluid and pore water are constructed by using a modified continuity equation in which a wave-generating source is included together with the modified momentum equations attached with fluid resistance terms, as shown in Eqs. (1)− (3).
where (x, z) is the Cartesian coordinate system, u and w are the velocity components in the directions of x and z, respectively, and γ v , γ x and γ z are the volume porosity and the projections of γ v in the x and z directions, respectively. t is time. g is the gravitational acceleration, and ρ is the density of fluid. p is the pressure, p=ρg z+p d , and p d is the dynamic pressure. τ is the viscous stress on the surface of the control volume. υ is the kinematic viscosity coefficient set as a constant of 0.010 1 cm 2 /s. β is the wave dissipation factor in the wave absorption zones. S * is the term of wave generation that is only defined at the source position, shown as Eq. (4), S is the flux density of the wave-generating source, and δ is Dirac's δ function.
Eqs. (5)−(7) are the resistance terms produced by porous materials, and they are defined as follows: M x and M z are the inertia forces, L x and L z are the linear resistances due to laminar flow, and T x and T z are the nonlinear resistances proportional to the square of velocity.
where C M , C L and C T are the empirical coefficients of inertia force, laminar resistance and turbulent resistance, respectively. D p is the median size of porous media in the seabed. Eqs.
(1)−(7) compose the unified mathematical model of the interaction between water waves and pore structures (including porous seabeds and permeable breakwaters). The different fluid domains in the model, i.e., the wave flow domain and the pore flow domain, are differentiated by the value of the porosity γ. When γ=1, the values of Eqs. (5)−(7) are all zero, and Eqs. (1)−(3) are restored as the well-known governing equations of the viscous incompressible fluid. When 0<γ<1, Eqs. (1)−(3) can model the pore flow, and γ is selected according to the actual porosity of the material in a permeable breakwater or a porous seabed. This model is relatively simple and avoids explicit boundary conditions at the interface between the external flow field and the pore flow. Moreover, the inertia force, laminar resistance and turbulent resistance generated by a permeable breakwater and seabed are adjusted by modifying the values of C M , C L and C T .

Boundary conditions
The governing equations mentioned above are constrained by the boundary conditions shown in Fig. 1. The boundary conditions of pressure and velocity should be satisfied on the free surface boundary, S f , of the fluid. The volume of fluid method (Hirt and Nichols, 1981) is used to track the free surface of the fluid control domain, where the change process of the fluid volume fraction function F over time is controlled by Eq. (8).
(8) Then, combining Eq. (1) and Eq. (8) into a conservative form leads to Eq. (9), (9) A donor-acceptor model is applied to calculate F. The pressure on the free surface is determined by linear interpolation or extrapolation. The flow velocity on the free surface is determined by the velocities between adjacent cells and the continuity equation.
The line source wave-generating method is employed to generate the desired incident waves. A source function wave-maker that can simultaneously produce two waves with opposite propagation directions is set up at a certain position in the calculation zone, as shown in Fig. 1. Added dissipation zones (Hinatsu, 1992) are arranged at both ends of the wave flume, where waves are damped by numerical dissipation effects due to the coarse grids. Open boundary conditions are applied to eliminate the wave secondary reflection. The fictitious damping forces based on Stokes's resistance law are taken into account in the z directions. Sommerfeld's conditions are adopted at the outflow boundaries so that the residual wave energy can spread to the exterior of the flume.
In addition to the abovementioned boundary conditions, there are impermeable solid boundaries denoted as S b in the model, such as the bottom of the wave flume and the impermeable seabed. A no-slip boundary condition, shown in Eq. (10), is adopted on the surfaces of these structures.
where s is the tangent direction of the solid wall boundary and n is the normal direction.

Numerical discrete schemes
The establishment of an efficient and stable numerical wave flume is the basis of numerical simulation for the interactions between waves and structures. Finite difference methods are applied to discretize the governing differential equations and the boundary conditions. The time derivatives are treated by the forward time difference scheme, the pressure gradient and stress gradient terms are treated by the central difference scheme, and the advection terms are treated by the combined scheme of central difference and the upwind method, with the weighting factors being 0.5. Staggered grids are adopted in the calculation area. The transient flow fields and pressure fields are obtained by the numerical algorithm of the SOLA-VOF method.
The determination of time interval Δt needs to satisfy the Courant−Friedrichs−Lewy (CFL) condition, as shown in Eq. (11), for computational stability.
where α is the weighting coefficient commonly taken as 0.5; |u| max and |w| max are the maximum absolute values of the flow velocity components in the x and z directions, respectively; and Δx i and Δz k are the mesh sizes. An initial value of 0.001 s is taken for the time step, which can be corrected by judging the CFL condition during the calculation progress.

Mesh convergence
Whether the discrete simulation results can truly reflect the physical phenomena described by the differential equation should be evaluated in terms of the verification and validation of the simulation results. In this study, Richardson's extrapolation method (Stern et al., 2001) is used to evaluate the numerical errors caused by the discretization of the computational domain. Considering the first-order precision extrapolation, three grid scales, Δx 1 =1 cm, Δx 2 =2 cm and Δx 3 =4 cm, are set in the computational domain with Δz all set to 1 cm, i.e., the grid scale ratio is 2. A third-order Stokes wave with incident wave height H=5.0 cm, period T=1.5 s, and water depth h=40 cm is selected for the working condition. Then, three solutions, S 1 , S 2 and S 3 , can be obtained, and the ratio is defined as R=ε 21 /ε 32 , which is used for mesh convergence judgment, where ε ij denotes the difference between S i and S j . Fig. 2 gives the numerical solutions calculated with the above three grids and the theoretical solutions for the thirdorder Stokes waves. This shows that the numerical solution can converge well to the theoretical solution when Δx 1 = 1 cm and Δx 2 =2 cm, and when Δx 3 =4 cm, the calculated results are relatively poor. The RMS errors of the three numerical solutions relative to the theoretical solutions are 6.64×10 −4 , 1.04×10 −3 and 1.93×10 −3 . The calculation yields 0<R<1, and it can be concluded that the calculated results converge monotonically to the grid scale. In the cases of three different grids, the spatial attenuation rates of wave height are 2.64%, 3.61% and 3.74% over the six wavelength ranges of the calculated region. Therefore, to ensure better accuracy, in the following numerical calculations, uniform rectangular grids with Δx=1 cm and Δz=1 cm are selected in the calculation area.

Wave making validation
To verify the stability of the wave generation of the numerical wave flume and the wave eliminating performance of the boundary, for the aforementioned wave elements, Fig. 3 presents the spatial distribution of wave profiles in different phases within one period in the calculation domain and the additional dissipation zone. In the calculation area (0.0<x/L<7.1, L=262.9 cm), the wave develops steadily in both the spatial and time dimensions, with the crest and the trough forming two horizontal envelope lines. In the wave elimination region on both sides (x/L<0.0, x/L>7.1), due to wave energy attenuation, the wave height decreases rapidly within two wavelength ranges and tends to zero at the open boundary. Good wave-making and wave elimination properties of the numerical wave flume can be observed.

Experimental set-up
In the physical model test, a seabed with limited thickness is designed as 500 cm long and 28 cm thick. The median grain size of the sand in the seabed, D 1 , is 0.033 8 cm. The porosity of the seabed γ is 0.4. The breakwater is arranged on the seabed with their midlines aligned to each other. The breakwater is 140 cm in bottom width and 30 cm in height, and its slope is 1:2 for both sides. The breakwater is composed of uniform rubble with a mean grain size D 2 =3.26 cm and a porosity γ=0.4. The wave height over the impermeable berm located at the front of the seabed is 5 cm, the wave period is 2.5 s, and the water depth over the mud line of the seabed is 40 cm. Several wave height gauges and pressure sensors are arranged for measuring wave surface elevation along the seabed and pore water pressure on the breakwater and seabed. Fig. 4 shows the layout of the structures and sensors. Fig. 5 shows a photo of the physical model test.
The geometric arrangement and the physical parameters in the numerical model are the same as those in the physical model test. A 3rd-order Stokes wave is numerically generated with the same factors as the experimental model. Due LI Yan-ting et al. China Ocean Eng., 2021, Vol. 35, No. 1, P. 26-35 29 to the larger particle size of the breakwater blocks, the nonlinear resistance plays a dominant role, and the influence of laminar flow resistance can be ignored. The resistance coefficients are taken as C M =1.67 and C T =16. In the seabed, inertial, linear and nonlinear drag forces are considered simultaneously with coefficients of C M =1.0, C L =0.12 and C T =2, respectively.

Resistance coefficients
The values of the resistance coefficients in the pore flow model are very important and directly affect the accuracy of the simulation result. The inertial force coefficient C M closely reflects the inertia of the fluid and the obstruction of the pore structure, i.e., the additional mass effect due to the change in velocity of the flow field within or around the structure. This can be expressed as C M =1+C m , where C m is the additional mass coefficient and C m ∈[0, γ/(1−γ)], then C M ∈[1, 1/(1−γ)] (Jung, 2004). For the selected porosity γ=0.4, C M ∈[1, 1.67]. As a result of the larger particle size of the blocks in the breakwater, the seepage flow velocity is also large, and the seepage field in the pore is unsteady. The effect of inertial resistance is very obvious; in this case, the inertia resistance coefficient can be taken as a larger value C M =1.67. In sandy seabeds, the influence of inertial force is relatively weak, and the lower limit value C M =1.0 can be used.
The effect of laminar flow resistance is evident at low Reynolds numbers; therefore, laminar flow resistance is only considered within the seabed in this model. Mizutani et al. (1996), in his study of wave-breakwater interactions, where the breakwater material was made of glass spheres with a diameter of 3 cm and the inertial, laminar and turbulent resistance was considered simultaneously, found that the dispersion degree of the distribution of C L grew larger with increasing KC number and that the relative importance of C L was slightly significant only if KC<10. Hur et al. (2008) took C L as 0.192 for a sandy seabed with a median particle size of 0.08 cm (note that as the expressions of the formulas in various studies are slightly different, the values of the resistance coefficient mentioned below are all converted to Eq. (6) or Eq. (7) in this paper). In this study, the particle size of bed sand is 0.034 cm; a smaller value is recommended, and C L =0.12 is finally determined.
C T is the nonlinear drag force coefficient, which reflects the viscous effect caused by fluid viscosity, whose value is affected by the fluid flow pattern, material particle size and surface roughness. Sakakiyama and Kajima (1992) determined C T by adjusting the reflection coefficient K r and transmission coefficient K t of the numerical model to make them approximately consistent with the value of the physical model when simulating the interaction between the waves and the permeable breakwater. Mizutani et al. (1996) found that C T was weakly dependent on the KC number, and C T was set as 0.23 in his study. Hur et al. (2008) assigned 60 and 25 as the values for a mound rubble breakwater (D=2.7 cm, γ=0.4) and sandy seabed (D=0.08 cm, γ=0.4), respectively. The value of the turbulence coefficient of resistance varies greatly from study to study. Referring to the above values and combining the results of multiple sets of tests, this study recommends taking 16 and 2 for the breakwater and seabed, respectively, so that the calculated results of the wave surface and the pressure inside the breakwater and seabed are in good agreement with the experimental results.

Numerical results validation
The interaction of waves, seabeds and permeable breakwaters has been simulated by using the present model. Comparisons are made between the numerical and experi-  mental results, and the numerical results are validated from three aspects, namely, the variation in wave height along the way, the pore water pressure on the breakwater, and the pore water pressure on the seabed. Fig. 6 shows a comparison of the wave time histories at several measuring points, and the wave heights and the profiles of the numerical and experimental results are in good agreement. Due to the integrated effects of wave reflection, shallow water and dissipation of the submerged breakwater, nonlinearity and asymmetry can be distinctly observed for the wave surface time history in front of and behind the breakwater. The trend of the average wave heights along the direction of propagation and the relative errors between the calculated values and the measured values are given in Fig. 7. Due to the reflection of the submerged breakwater and the permeation of the seabed, the wave height at measuring point No. 2 is significantly reduced, and it begins to increase due to the effect of shallow water above the wavefacing slope of the breakwater; then, it decreases at the top and behind the breakwater with the dissipation of wave energy. The relative errors of the wave heights are large, reaching 14.98% (No. 3) and 13.71% (No. 5), because of the complicated nonlinear evolution of the wave surface over the wave-facing slope and the crown of the breakwater.   Fig. 8 depicts the time-history curves of the wave-induced dynamic pore water pressure at each measuring point on the symmetry axis of the breakwater numerical simulation and the experimental test. The corresponding mean values of oscillation amplitude (denoted by P A ) and relative errors are listed in Table 1. Due to the dissipation of wave energy by breakwater materials, the amplitude of the pore water pressure in the submersion breakwater decreases gradu- Fig. 8. Time history of dynamic pore water pressure at the measuring points on the breakwater. ally with depth. Fig. 9 shows the wave-induced dynamic pore water pressure progress at four measuring points along the depth at the midline of the seabed. Table 2 lists the calculated and test values of the oscillation amplitude and their relative errors. In Fig. 10, the calculated and test values of the average pore water pressure amplitude at each measuring point in the seabed are compared, and the positions of the sampling points correspond with those in Fig. 4 one by one. The calculated value at point No. 19, directly below the submerged breakwater near the seabed surface, is smaller than the actual measured value. The calculated values at most of the remaining points match the experimental values well. The amplitude of pore water pressure in the seabed decreases gradually with depth and presents an obvious nonlinear distribution. The gradient of fluctuating pore water pressure near the seabed surface is obviously larger than that in the deeper part of the seabed in the vertical direction, and the amplitude at the deeper part of the seabed essentially decreases linearly. Therefore, the soil particles in the surface layer of the seabed under the breakwater are prone to instability due to wave action.

Interstitial flow fields
The seepage flow in the breakwater and seabed induced by waves has an impact on the local stability of the permeable breakwater and may change the geometrical and mechanical properties of the seabed as well, thus threatening the structural safety. The correctness of the numerical results has been verified in Section 4.3, and a more detailed analysis is carried out below. It should be noted that the scale of the permeable seabed, in other words, the solid-wall boundary (horizontal bottom and impermeable seabed on both sides) existing in the model, presents a certain degree of influence on the seepage at the impermeable boundary and has less effect on the seepage near the breakwater, the wave movement outside the structure and the dynamic pore water pressure in the seabed and breakwater, which can also be seen in the subsequent Figs. 11 and 12.
The instantaneous flow fields around the seabed are shown in Fig. 11, where the individual arrows denote the velocity vectors and the background color identifies the velocity magnitude. The scale of the pore flow velocity in the seabed is far less than the flow velocity in the wave field. From the pattern of streamlines, an obvious phase difference is observed between the wave flow and the seepage in the seabed. The same phenomenon was found by Hur et al. (2010) in his study, who suggested that the phase difference was caused by the laminar flow resistance in the seabed arising from the "pumping" effect (Elliott and Brooks, 1997), the movement of fluid through the seabed generated by steady periodic dynamic pressure over the seabed surface. In Fig. 11a, a wave crest is moving towards the breakwater. The surface elevation is relatively low above the seaward slope. The fluid overflows the breakwater from the leeward to the seaward, and a vortex forms near the seaward toe of the breakwater. The seepage inside breakwater flows from the leeward slope to the seaward slope. The direction of the seepage in the seabed is from leeward to seaward. Fig. 11b shows a state in which the wave crest reaches above the toe of the breakwater. The water flows up along both the seaward and the leeward slopes. Due to the shallow water effect, the wave height and the flow velocity above the seaward slope both increase. Part of the seepage in the breakwater flows from the leeward slope to the crest, and the other part flows from the seabed to the breakwater and converges on the seaward slope. The seepage in the seabed flows from the mud surface front and behind the breakwater to the breakwater bottom and then flows into the breakwater. In Fig. 11c, the crest reaches the top of the breakwater, and the wave height and the flow velocity increase. Strong shear flow forms on the seaward slope of the breakwater. The seepage in breakwater flows from the seaward slope to the leeward slope. The seepage in the seabed flows from the bed surface in front of the breakwater to behind the breakwater, and a vortex forms near the mud surface in front of the breakwater. In Fig. 11d, the wave crest  LI Yan-ting et al. China Ocean Eng., 2021, Vol. 35, No. 1, P. 26-35 runs over the top of the breakwater, with the wave energy decaying and the wave height decreasing. The maximum flow velocity occurs at the top. Part of the seepage flows from the seaward slope, and the other part infiltrates into the seabed. The pore water in the seabed flows from the vicinity of the toe to the front and rear of the breakwater. The vortex in the seabed behind the breakwater is formed by the blocking of the impermeable seabed. In Fig. 11e, the wave damped by the permeable breakwater and seabed spreads to the rear area and forms a strong downward shear flow on the leeward slope of the breakwater, and a small vortex forms at the rear toe of the breakwater. Part of the water flow in the breakwater seeps into the seabed from the leeward slope, and the other part seeps from the leeward slope to the seaward. The interstitial flow in the seabed oozes from the underneath of the breakwater to the mud surface in front of and behind the breakwater. Then, the wave moves away, and the next cycle begins in Fig. 11f.
Throughout the entire cycle of the wave passing through the submerged breakwater, the maximum velocity occurs over the upper part of the seaward slope, the crown and the shoreward top corner of the breakwater; thus, these parts will suffer from greater wave loads. On the two slope faces of the submerged breakwater and both sides of the seabed surface near the breakwater, stronger shear flow occurs with larger velocity gradients. Vortexes tend to form at the foot of the breakwater and near the bed surface in front of the breakwater, and protection should also be strengthened at these positions. Fig. 12 displays the instantaneous distributions of the wave-induced dynamic pressure at different phases that correspond to the phases in Fig. 11, where the background color denotes the magnitude of pressure. The pattern of the dynamic pressure distribution shown in Fig. 12 is consistent with the seepage in Fig. 11, i.e., the flow moves from the high-pressure area to the low-pressure area. As the wave propagates to the positions in Figs. 12b and 12e, the dynamic pressure in front of the breakwater is reduced compared with that above the impermeable bed surface due to the permeability of the seabed. As the wave continues to travel backwards, the pressure in the external flow field starts to increase due to the presence of the breakwater. Since the wave energy has been partially dissipated, the wave height decreases and the wave pressure decreases accordingly behind the breakwater. In addition, the pressure gradient on the surface of the pore structure and the interface between the seabed and breakwater is significantly larger, and therefore, the seepage flow is stronger here. The pressure distribution also reveals a phase difference between the external wave field and the internal flow of the pore structure. Fig. 13 displays the spatial distribution of the maximum dynamic pore water pressure normalized by ρgH inside the permeable breakwater and the seabed. The maximum dynamic pore water pressure in the permeable breakwater is distributed in the middle and upper parts of the seaward slope and the crest, with a value of 0.59. Waves mainly break in these positions, and the wave energy is dissipated in large quantities. The pore water pressure near the back face of the breakwater is also relatively high. The pore water pressure in the breakwater basically decreases along the inner normal direction of the two slope faces. In general, the maximum dynamic pore water pressure in the breakwater near the seaward side is greater than that of the other side. Where the distribution value of the pore water pressure gradient is large, stacking stones are prone to instability. Fig. 14 illustrates the distribution of the maximum fluctuating pore water pressure in the submerged breakwater under the same working conditions when a non-impermeable seabed is considered. By comparing Fig. 13 with Fig. 14, it can be seen that the permeability of the seabed has a significant impact on the pressure distribution within the breakwater. In both permeable and impermeable seabeds, the value of the maximum dynamic pore water pressure in the breakwater of the latter is generally greater than that of the former, and the pressure distribution of the latter basically varies along the horizontal direction, with a small pressure change rate along the vertical direction, which is much smaller than that of the former.  The maximum dynamic pore water pressure in the seabed gradually decreases along the depth as a whole and shows a decreasing trend in the wave direction under the breakwater. The maximum dynamic pore water pressure in the seabed is located under the seaward toe of the permeable breakwater, with a value of approximately 0.35. The dynamic water pressure gradient near the seabed surface at the bottom of the breakwater is large, which is likely to cause seepage between the seabed and the breakwater and easily cause the loss of soil particles on the seabed, adversely affecting the stability of the foundation and threatening the safety of the breakwater.

Conclusions
In this study, the interactions among waves, permeable breakwaters and seabeds are calculated by applying a twodimensional difference numerical wave flume based on an improved N−S equation. The calculated results are in good agreement with the experimental values. In addition, the results are analyzed in detail. The main conclusions are summarized as follows: (1) The calculated and tested results of the wave height in front of and behind the breakwater and the pore water pressure inside the permeable breakwater and the seabed are compared, confirming the reliability of the calculated results. The model can be used for a numerical simulation of the wave−permeable breakwater−seabed interaction.
(2) Analyses of the seepage flow fields and the pressure distributions in the breakwater and the seabed show that strong shear flows and relatively high pore water pressure gradients exist on the slope and the crown of the breakwater, and these parts are prone to structural instability.
(3) The gradient of the oscillation scale of the dynamic pore water pressure in the seabed is nonlinearly distributed in the vertical direction. The pore water pressure gradient pressure close to the seabed surface beneath the breakwater bottom is relatively high, especially that adjacent to the toe of the breakwater at the seaward side, which is prone to cause strong seepage.
(4) The permeability of the seabed has an obvious influence on the spatial distribution of dynamic pressure on a permeable breakwater. When the seabed is impermeable, the spatial distribution of hydrodynamic pressure mainly changes along the horizontal direction. When the seabed is permeable, the vertical change rate of the dynamic water pressure increases significantly, which is unfavorable to the stability of the seabed and the breakwater.
At present, the determination of resistance coefficient values varies considerably in similar studies, but its importance is self-evident. Subsequently, the variation law of experimental factors on the selection of the resistance coefficient will be studied in detail, including pore material characteristics, breakwater types, working conditions and model scales, to further develop the engineering applications of such a model.