Investigation on A Resistance-Type Porosity Model and the Experimental Coefficients

The mathematical model based on the Volume-Averaged/Reynolds-Averaged Navier—Stokes (VARANS) equations has been adopted in recent years to generally simulate the interaction between waves and porous structures. However, it is still hard to determine the two experimental coefficients (α and β) included in VARANS equations. In the present study, VARANS equation is adopted to describe the flow inside and outside the porous structures uniformly, with applying the volume averaged k − ε model to simulate the turbulence effect. A new calibration method is used to evaluate the accuracy of numerical simulation of wave motion on porous seabed with different coefficients, by taking the wave damping rate as an index. The calibration is achieved by completing a simulation matrix on two calibration cases. A region can be found in the parameter space to produce a lower error, which means that the simulation results are better consistent with the experimental results.


Introduction
The simulation of the flow in porous structures by numerical methods and computational resources originated in the late 1900s. Many researchers have studied the establishment and numerical methods for coupling free and porous media flows. Van Gent (1995) added a resistance source term into the Navier-Stokes equation describing the effect of the porous media. This source resistance term can be expressed as the extended Darcy-Forchheimer equation which contains an inertia term, a linear resistance term and a quadratic resistance term. These three terms exert forces on the fluid due to acceleration, drag and friction. The resistancetype porosity models developed by Van Grent (1995) solved some of the wave-porous structures interaction problems. However, the limitation of this model is that it does not consider the turbulent effect of the flow field inside and outside the porous media. Liu et al. (1999) adopted the formulation in van Grent (1995) by further considering the turbulence outside the porous structures. The turbulent boundary layer adjacent to the porous wall is modified by including the effect of percolation velocity along the porous boundary. Hsu et al. (2002) developed a mathematical model based on the Volume-Averaged/Reynolds Averaged Navier-Stokes (VARANS) equations to describe the flow both inside and outside the porous medium. In order to obtain the VARANS equations, intrinsic volume average is applied to the RANS equations, and some additional terms including residual stress term and two facial integral terms occurred. Hsu et al. (2002) modeled these terms as extended Darcy-Forchheimer relationship collectively which can be expressed as: These two resistance coefficients (a, b) are crucial in modeling the interaction between simulated waves and porous structures, which must be correctly determined. Engelund (1953) formulated the relation between these resistance coefficients and the porosity, the viscosity and grain diameter for steady-state flow, as Burcharth and Andersen (1995) and del Jesus et al. (2012) did in their work. Van Gent (1995) formulated similar expressions which added the effect of oscillatory flows to the expressions in terms of KC number. The latter formulations were applied in the model presented by Liu et al. (1999).

α β
More troubling, these two resistance coefficients are also related to the well-known experimental coefficients and . It should be noted that different values were given by researchers depending on how the coefficients were determ-α β α β α β ined or calibrated, how the resistance coefficients a and b were formulated and even how the turbulence effect was dealt with by the mathematical model. The formulation by Engelund (1953) yielded values of and differ from those obtained in the formulation by van Gent (1995). Jensen et al. (2014) proposed a new method to determine and . They applied three cases of dam break through porous media with different Reynolds numbers. By comparing the free surface between numerical and experimental results, and were calibrated to a reasonable range. However, their error estimation process was too complicated because the water surface location is considered not only along the horizontal distance but also along the time steps.
In this study, VARANS equation was adopted to describe the flow inside and outside the porous structures uniformly, with applying the volume averaged model to simulate the turbulence effect. In addition, a new calibration method was applied that the wave damping rate can be regarded as an index to evaluate the accuracy of numerical simulation of wave damping over porous seabed when using different and values. Fig. 1 showed the definition of porous media. The surface S, which is used to create the averaging volume V, may include both solid phase and fluid phase. S is defined by a circle with the radius r 0 . The actual volume of the fluid phase V f may vary with the porous media based on the position of the averaging volume while volume V is constant. The macroscopic length scale L and the pore length scale l are also shown in this figure.

⟨⟩ ⟨⟩ f
The volume averaging procedure was applied to transform the RANS equations to the VARANS equations and the length scale constraint was given as l<< r 0 <<L. The superficial volume average denoted by was defined as the average over the entire volume. The intrinsic average denoted by was defined as the average over the fluid volume. The superficial and intrinsic volume averages of a scalar, vector, or tensor denoted by B were defined as: (2) The relation between these two averages was: where n was the porosity given by: In order to obtain the VARANS equation, the Navier-Stokes equations were ensemble averaged to obtain the Reynolds Averaged Navier-Stokes equations (RANS) by applying the time-average operator. Then, the volume-average operator was applied to the RANS equations to obtain the Volume Averaged/Reynolds Averaged Navier-Stokes equations (VARANS): ∂ where f i along with the inertia term was the extended Darcy-Forchheimer resistance term defined as: α β The first term on the right hand side of Eq. (8) was the linear resistance term, and the second term was the quadratic resistance term. and were experimental coefficients to be determined in advance, which was the main concern of this paper. The last term was the inertia term. c m =0.34(1-n)/n was the added mass coefficient which accounts for the transient interaction between grains and water. It was obvious that n=1, c m =0 were outside the porous media, then VARANS equation automatically changes back to RANS equation.
In Eq. (7), was the volume-averaged eddy viscosity, the balance equations of and can be obtained by taking the volume averaging over the standard k and equations which turn to: . and are the additional turbulence sources due to the presence of porous materials. They represented the small-scale turbulence with a scale smaller than the average volume. Nakayama and Kuwahara (1999) proposed closure schemes for small-scale turbulence based on the numerical simulations of the flow through a series of square rods: 3 Calibration method α β These two experimental coefficients and in the extended Darcy-Forchheimer equation existed only in the resistance-type porosity models without any physical explanation. Although the recommended values are given in the references, other researchers can not directly use these values without calibration if different formulas of resistance coefficients and different methods of dealing with turbulence effects are taken into account.
When the wave propagated over a porous seabed, the wave height decreased due to the flow resistance inside the porous bed. Assuming that the free surface varied with distance and time as follows: σ where a was the wave amplitude and was the angular frequency. The wave number k w , a complex variable, can be written as k w =k r +k i , where k r represented the change of wavelength and k i represented the damping of wave height. Considering that the value of k r was small, the relation between wave height at the location '1' and '0' with a distance of x can be described as: When a resistance-type porosity model is used to simulate α β α β the wave motion on the porous seabed, different selections of and values can lead to different results on the wave damping rate k i . With reference to Hsu et al. (2002), Lara et al. (2010) anddel Jesus et al. (2012), the calibration was achieved by completing a simulation matrix on two calibration cases which come from the experiments by Savage (1953) and by Sawaragi and Deguchi (1992). The two coefficients were varied as = [100,300,800,1200,2000] and =[0.5, 1.0, 1.5, 2.0, 2.5, 3.0]. The error between simulation and experiment was estimated according to the difference between the simulated wave damping rate k i(c) and the measured wave damping rate k i(m) : A region can be found in the parameter space to produce a lower error, which meant that the simulation results were better consistent with the experimental results.

Sawaragi and Deguchi's experiment
The experiment of Sawaragi and Deguchi (1992) was carried out in a two-dimensional wave tank with length of 30 m, width of 0.7 m and depth of 0.9 m. The permeable bed was made at the horizontal bottom of the wave tank. The length and the thickness of the permeable bed were 3.5 m and 0.15 m, respectively. Rubble stones with an average diameter of 3.07 cm were used to construct the permeable bed. Since the porosity was not provided in their paper, we use n=0.4 estimated by Chang (2004). The water depth above the permeable bed was 0.15 m. The wave period T=1.0 s and the wave height H i =0.0358 m. Fig. 2 shows the sketch definition of the numerical model according to the experiment. The permeable bed was located 2 m away from the wave maker. A total of 61 m of wave height was evenly arranged above the permeable bed. A close distance between gauges of 0.1 m was set to ensure small error of the curve fitting wave height damping. Fig. 3 and Fig. 4 show the wave variations above the impermeable and permeable seabed respectively. It can be seen from Fig. 3 that the wave height hardly changed with the distance when the seabed was impermeable. It meant that the wave damping was negligible. On the other hand, Fig. 4 shows that the wave height decayed obviously when the seabed was per- Fig. 2. Performance of the numerical wave tank on wave motions over permeable seabed. meable, which meant that a significant wave damping was introduced. α β Sawaragi and Deguchi (1992) gave wave heights over distance and the wave damping rate was about 0.1463. By completing a simulation matrix, where the two experimental coefficients were varied as = [100,300,800,1200,2000] and = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0], a total of 30 wave damping rates were computed as shown in Table 1.
As mentioned above, the wave damping rate was affected by the experimental coefficients and . Taking part of the results as an example, namely, =300 and varied as = [0.5, 1.0, 1.5, 2.0, 2.5, 3.0]. Fig. 5 showed that the larger gave rise to more obvious the wave height decay. Fig. 6 showed that led to a similar effect on wave damping with . It turned out that, in Sawaragi and Deguchi's case, wave damping rate was positively correlated to both the two experimental coefficients and . The error between the simulated and measured wave damping rate k i(c) was presented in Fig. 7 as the contours of the parameter space. It was found that not one unique combination of and would give a minimum error. On the contrary, a slanted band with small error is found in the parameter space, and the error increases gradually on both sides of the band. The contours are aligned along the axis representing both and , which meant that the two experimental coefficients affected the wave damping rate in the same way.    The brighter color regions in the contour represent a lower error parameter space, but it was still too large to determine a meaningful combination of the two coefficients. Therefore, the Savage experiment was analyzed as another calibration case to further limit the range of and .

Savage's experiment
The experiments were performed in a wave tank with the length of 29.3 m, the width of 0.46 m and the depth of 0.61 m. The permeable bed was made at the horizontal bottom of the wave tank. The length and the thickness of the permeable bed were 18.3 m and 0.3 m, respectively. The average diameter of porous media was 3.82 mm. Since the porosity was not provided in Savage's paper, the n=0.3 was used in this work, which was estimated by Liu and Dalrymple (1984), Gu and Wang (1991) for the same experiment. The water depth above the permeable bed was 0.229 m. The wave period was T=1.27 s and the wave height was H i =0.054 m. It can be seen that the average diameter of the porous media used by Savage was much smaller than that of Sawaragi and Deguchi (1992), and the porosity decreased in the same way. The flow in the porous media slowed down and the turbulence effect tended to be weak. Theoretically, when the wave height attenuation was simulated by a resistance-type porosity model, the linear drag term played a more important role than the nonlinear drag term.  Table 2.
The error between the simulated and measured wave damping rate k i(c) was shown in Fig. 8 as contours over the β α β parameter space. Different from Fig. 7, most of the outline in Fig. 8 was arranged along the direction of the -axis, indicating that the error mainly depended on the value of , the effect of can be therefore neglected.
Because the brightest color in the contour, i.e. the lowerror space ( <300), was too close to the left boundary, two additional values of ( =50 and =200) were added to the simulation matrix as shown in Table 3. Fig. 9 shows the refined contour with ranging from =50 to =300. The results showed that when was larger than 100 and smaller than 200, the vertical band provided low errors and the error increased significantly when was smaller than 100.

α α
It was noteworthy that another difference was found between Sawaragi and Deguchi's case and Savage's case. In Savage's case, it can be seen that the wave damping did not increase along . On the contrary, it decreased when varied from 200 to 2000. The interaction between waves and porous seabed seemed complex. The wave damping rate Fig. 8. Contours of the error between simulated and experimental wave damping rate for Savages case. Fig. 9. Contours of the error between simulated and experimental wave damping rate for Savages case (refined part). was not always positively correlated to these two experimental coefficients, because the percolation in porous structures was weakened when these two coefficients increased. On the other hand, it indicated that these two coefficients were kind of constants and meaningless out of the resistance-type porosity models. According to the calibration process presented in Sections 3.1 and 3.2, it was found that a common area with varied as [100,200] and varied as [1.6, 2.2] provided an optimal solution for the two calibration cases. This solution, especially the value, was smaller with reference to Jensen et al. (2014) which was most likely because there was no turbulence model adopted in his work and thus the turbulence effect was described via the Darcy-Forchheimer resistance term. Particularly, =100 and =2 were considered to be the optimal combination of the coefficients and would be used in the following sections.

Validation cases
4.1 Wave height over distance above porous seabed α β By applying optimal coefficients =100 and =2 in the porosity model, the wave height along the distance corresponding to the experiments of Sawaragi and Deguchi (1992) and Savage was numerically simulated. Fig. 10 and Fig. 11 showed the comparisons between numerical results and experimental results. It can be seen from these figures that both numerical results were in good agreement with the experimental results, even considering that the permeable seabed of Savage was much longer than that of Sawaragi and Deguchi (1992).

Wave patterns in front of a perforated caisson on a rubble bed
The wave pattern in front of a perforated caisson on a rubble bed was experimentally studied, and the wave reflection was also studied. The depth of water was 0.4 m. In the middle of a rubble-mound foundation, caissons with perforation ratio of 40%, wave chamber length of 0.2 m and 0.3 m were applied. The foundation was 0.2 m high and the top of the foundation had a length of 95 cm. Two slopes of 1:2 were set on the front and back side of the foundation respectively. The average diameter of constructed stones was 1.1 cm and the porosity was 0.3. In order to collect free surface data and calculate the wave reflection coefficient, five gauges with different spacing were installed in front of the caisson. According to the experiment, the numerical wave tank was performed as Fig. 12.
When the wave height was h i =0.06 m, wave period is T=1.4 s and the length of wave chamber is b c =0.3 m, the wave reflection coefficient calculated from experimental data was about 0.37, and which calculated by numerical simulation was 0.42. Fig. 13 showed the time series of water surface at the locations of wave height gauges #1-#5 obtained by numerical simulation and experiment. It was illustrated in this figure that the numerical results were in good agreement with experimental results of both wave periods and wave lengths at all other locations, except for some deviations at gauge #3. Fig. 14 showed the comparison of the time series of water surface when the wave height h i =0.08 m, wave period T=1.2 s and the length of wave chamber   ZHAO Pei-hong et al. China Ocean Eng., 2019, Vol. 33, No. 4, P. 468-476 b c =0.2 m. The wave reflection coefficient calculated from the experimental data was about 0.44, and that calculated by the numerical simulation was 0.43.

Run-up and run-down pressure in rubble-mound foundation and backfill material
The rubble-mound foundation set beneath the caisson not only affected the wave pattern and the wave reflection coefficient, but also affected the wave force acting on the caisson structure both in horizontal and vertical directions. Especially when backfill permeable material exists in the wave chamber, run-up and run-down pressure in porous me-dia play an important role in vertical forces. The experimental study on the vertical wave force on a perforated caisson setting on a rubble-mound foundation was carried out by Wan et al. (2010). The experimental set-up was almost the same as mentioned in Section 4.2, except that the height of rubble-mound foundation was reduced to 0.10 m. As shown in Fig. 15. The arrangements of pressure gauges were shown in Fig. 15 and Fig. 16. Six gauges were arranged on the upper side of the chamber and five gauges were arranged on the lower side, which are expressed by circles and triangles respectively. When the width of wave chamber was reduced from 0.30 m to 0.20 m, some of the upper gauges would become redundant and no data would be collected. On the other hand, the pressure acting on the lower side of the chamber floor was negative and reduced with the distance from the left boundary. Fig. 18 shows the negative pressure amplitudes. It can be seen from this figure that the pressure acting on the upper side of the chamber floor was negative and increased with the distance from left boundary. The pressure acting on the lower side of the chamber floor was positive and reduced with the distance from left boundary. After all, the change of the run-up and run-down pressure near the left boundary of the chamber floor were more significant than those far from the left boundary. In addition, the numerical results were in good agreement with the experimental results and showed similar variation trends. Fig. 19 and Fig. 20 show the positive and    negative pressure amplitudes when h i =0.08 m, T=1.2 s and b c =0.2 m. It was demonstrated that the variation of both numerical and experimental results were basically consistent, and the simulated data was in good agreement with the experimental data.

Conclusions
The mathematical model based on the Volume-Averaged/Reynolds-Averaged Navier-Stokes (VARANS) equations was used to simulate the interaction between waves and porous structures. However, it was difficult to determine the two experimental coefficients included in VARANS equations. In this paper, a new calibration method was ap-α β plied to evaluate the accuracy of numerical simulation of wave damping over porous seabed under different values of and values, by taking the wave damping rate as an index. According to experiments of Savage (1953) and Sawaragi and Deguchi (1992), two calibration cases were used to complete the simulation matrix and accomplish the calibration. Then, several validation cases were carried out. The results show that the wave patterns in front of a porous structure and the wave-induced run-up and run-down pressure in porous materials were in good agreement with the experimental data. Generally, the whole calibration and verification process was clear and easy to be carried out. The new calibration method can be used for reference.