Porosity Effects on the Dispersion Relation of Water Waves through Dense Array of Vertical Cylinders

: There is growing interest for water ‐ wave flows through arrangements of cylinders with application to the performance of porous marine structures and environmental flows in coastal vegetation. For specific few cases experimental data are available in the literature concerning the modification of the dispersion equation for waves through a dense array of vertical cylinders. This paper presents a numerical study of the porosity effects on the dispersion relation of water waves through such configurations. To this aim, the sloshing problem in a tank full of vertical cylinders intersecting the free surface is studied using the finite element method, and the influence of the porosity on the wave number is quantified. On the basis of numerical results, a new modification of a dispersion relation for porous medium is suggested based on a wide range of collected data. Moreover, the domain of validity of this new dispersion relation is examined considering the number of cylinders and the extrapolation to the infinite medium.


Introduction
For many years, the scientific community dedicated significant efforts to the understanding and modeling of the wave-structure interaction; see, e.g., Mei [1] and Molin [2]. Among the structures, many studies are focused on cylinders, since they constitute a basic component of rigs and structures employed in ocean and coastal engineering applications. Many experimental works (e.g., Goda [3]), mathematical models and numerical methods (e.g., John [4], Yang and Ertekin [5]) are developed to understand the reference problem of water waves interacting with a single, bottom mounted or floating, surface piercing, vertical cylinder; see also Young [6], Sabunku and Calisal [7]. A review can be found in McNatt [8].
However, with the development of complex ocean structures, such as very large floating structures (Wang and Tay [9]), arrays of devices for the outtake of marine renewable energies (Falnes [10], Balitsky et al. [11]), and for porous coastal and harbor defenses (Arnaud [12]), the problem has been expanded to more complex geometrical configurations, where dense arrangements of vertical cylinders are involved. Similarly, there is growing interest concerning similar problems in environmental flows in vegetation, such as the propagation of water waves in mangrove forests (e.g., Massel [13], Mei et al. [14]).
The above configurations, of course, lead to a large increase of the difficulty, which has been addressed by using various strategies. An important body of the literature is dedicated to the study of water waves propagating through dense arrays of cylinders. If it is not possible to be exhaustive, one should emphasize the studies of Linton and Evans [15], Maniar and Newman [16], McIver [17] or Kagemoto et al. [18]; see also Belibassakis et al. [19].
From the point of view of numerical methodology, the direct calculation of water waves through dense arrays of cylinders in complex configurations is computationally demanding. Of course, recent improvements both in computational resources and in numerical optimization techniques offer possibilities to solve this problem when the degree of complexity remains reasonable (Fu et al. [20]). Thus, another possible approach is to rely on extensions of the potential flow theory, whereby the cylinder domain is treated as a homogeneous porous medium. The dispersion equation, though, is modified to take the extra inertial and dissipative terms into account. Such formulations, for instance, are suggested by Sollitt and Cross [21], Madsen [22] or Yu and Chwang [23]. These modifications of the dispersion were used recently to model the reflection, transmission (Arnaud et al. [24]) and diffraction (Arnaud et al. [25]) of water through an idealized porous medium of various porosity and specific surface.
In the meantime, by implementing an experimental approach relying on the sloshing response of a tank filled with homogeneous vertical, surface piercing cylinders, Molin et al. [26] questioned the formulation of the dispersion equation to be considered and introduced a new expression. An illustration of the tank used in the above work, which filled with a regular pattern of vertical cylinders, is presented in Figure 1. Although the sloshing problem finds applications to engineering mechanics, the eigenvalue problem of fluid motion within tanks containing arrays of surface piercing obstacles provides also useful information for environmental maritime problems. The present work is a numerical extension of the study by Molin et al. [26]. The resonant sloshing in tanks with internal (sloshing-suppressing) structures has been discussed by various authors; see, e.g., Faltinsen and Timohka [27]. Although the solution of the nonlinear problem is very important still the solution of the linearized problem provides useful information concerning the natural sloshing modes and frequencies for generalizing the method (Timohka [28]). Since the structures (cylinders and tank walls) contain vertical boundaries the 3D problem can be reduced to 2D governed by the Helmholtz equation on the horizontal plane. An important novelty of the present approach is that it is an easily applicable method to more general shapes of tanks, with arbitrary distributions and general cross section characteristics of the structures in the tank.
A numerical solver, providing the eigenmodes of the tank, based on a finite element approach, is developed, allowing one to extend the analysis of the influence of the porosity τ of the medium (defined as the ratio of water volume in terms of the total volume of the domain) on the dispersion equation.
The paper is organized as follows. In Section 2, the equations of the problem and the FEM numerical procedure for the solution are presented. A comparison of the tool with experiments by Molin et al. [26] is also performed here, to provide an insight of the performance of the present model. Next, Section 3 focuses on the analysis of the effect of porosity on the equivalent dispersion equation, or equivalently to the equivalent bulk wavenumber of water oscillations developed in the tank. A modification of this equation is introduced. Finally, results are summarized in Section 4, where also conclusive remarks are provided.

Mathematical Model for the Dispersion Relation
The sloshing problem of water in a tank filled with water up to a level and containing the dense arrangement of vertical cylinders is considered using potential flow theory, in conjunction with the assumption of small amplitude waves. Since the flow is irrotational, the velocity is obtained as the gradient of a 3D velocity potential Φ: (1) Equation (1) in conjunction with the incompressibility of the water, leads to the Laplace equation: We consider the condition of impermeability at the bottom of the tank (z = −h): and the linearized free surface boundary condition at z = 0: In the case of harmonic motion at given frequency ω, exploiting the fact that the walls of cylinders and tank boundaries are vertical the evanescent modes are not excited and the wave field Φ is represented by: where is the vertical profile of the potential field and the corresponding 2D velocity potential on the x-y plane. With this form of potential, Equation (2) leads to: where a prime denotes differentiation with respect to the involved variable and Δ stands for the twodimensional Laplace operator. This means that that the quantity is constant everywhere in the water domain. The latter is denoted by , and thus, represents the wave number on the horizontal plane. Consequently, Equation (6) is written: Furthermore, from Equation (7) we obtain that the 2D potential satifies the Helmholtz equation on the horizontal plane: valid in the water domain. The solution of Equation (7) satisfying the boundary condition (3) is: which substituted in the potential Φ x, y, z, t , and used in the free surface condition (4) leads to: In the sequel is considered constant everywhere in the domain, defined as .
It is noticed here that the associated wavelength , does not correspond to the water wavelength , which for the first sloshing mode is twice the tank length L: as implied by the fulfillment of the side boundary conditions on the vertical tank walls. The objective of this study was to determine a relation between and , in order to associate the water wavelength to its frequency in the porous medium In Section 2.2 a numerical sloshing tank was modeled by using FEM to establish the above relation.
A first set of calculation, presented in Sections 3.1 and 3.2, will be useful to optimize the cylinders arrangement of the numerical sloshing tank in order to cover the widest range of cases and reduce computation time. Then, many cases were tested covering a porosity of τ = 0.3 up to 0.95 and a novel formulation of was proposed. The results were also compared against the modification of the dispersion relation proposed by Molin et al. [26], which was also expressed by Equation (15) with G(τ) = τ 1/2 .

FEM Numerical Scheme
The present numerical method was based on the triangulation of the horizontal domain on the free surface and application of standard first-order FEM to the solution of the eigenvalue problem. As an example of verification, we first considered the case of sloshing modes in a circular container with vertical cylinders (poles) for which analytical solutions were available; see e.g., Timokha [28].
A particular example is presented in Figure 2 for a container of unit radius with a single eccentric circular pole of radius 0.5 m centered at x = 0.5 m, considered by Timohka ( [28], Section 4 and Figure  3). A sketch of the triangulation corresponding to a mesh of 4693 elements, which was proved enough for numerical convergence, is shown in Figure 2a. The lowest eigenfunctions corresponding to the symmetric and antisymmetric modes are plotted in Figure 2b,c, respectively and were found in very good agreement with the above analytical solution. Finally, the result of calculation by the present FEM concerning the lowest eigenvalue corresponding to the symmetric mode in the same container for various values of the radius Rin of the inner circular pole is illustrated in Figure 3. The latter quantity is reported to be not well predicted by the asymptotic formula provided by Faltinsen and Timohka ( [27], Equation (4.185)). However, the present results are found in good agreement with analytical predictions obtained by the variational Trefftz method by Timohka ( [28], Figure 2).

Algorithm
The present study requires FEM calculations using a large amount of meshes. Consequently, it was necessary to automatize the mesh generation of a numerical sloshing tank. The mesh generator was implemented based on the Delaunay triangulation. Among other information, the algorithm was based on the following more important input data:


The tank size length and width (L and W);  The number of cylinders (therefore Nx and Ny for the number in both direction);  The porosity of the medium;  The arrangement of the cylinders;  The level of refinement.
Then, the algorithm follows the following steps: 1. Calculate the cylinder radius accounting for the target value porosity (if necessary); 2. Calculate the position of the center of cylinders accounting for their number and the target arrangement; 3. Generate an initial mesh; 4. Optimize the triangle quality by slightly moving nodes; 5. Refine the mesh.
Steps 4 and 5 are repeated as many times as chosen by the user with the input "level of refinement".

Arrangements
Four different arrangements are automatized as shown in Figure 4. The experimental arrangement Figure 4a corresponds to the tank used for the experimental campaign conducted in Marseille and described in Molin et al. [26]. For aligned and misaligned arrangements Figure 4b,c, the distance between two successive cylinders remains the same lengthwise and widthwise (not necessarily the same distance in both directions). Furthermore, the distance to a wall is half the distance between two cylinders. Finally, the random arrangement Figure 4d follows no rules except the cylinders cannot intersect each other.

Repeatability of the Pattern
If several aligned or misaligned meshes are repeated lengthwise or widthwise, the arrangement remains the same as shown in Figure 5. In particular, results obtained for a given mesh can be extrapolated to an infinite medium.

Numerical Solver
The numerical solver was based on the finite element method (FEM) to solve the Helmholtz Equation (8); see e.g., Johnson [29]. Let us consider the weak form of the Helmholtz Equation (8) Δ , where , are the linear pyramid functions equal to 1 for node number i, and 0 elsewhere and Ω is the water domain on the horizontal plane. Denote the approximation of the field as where the nodal values in the form of a global vector are ⋮ .
Using the FEM approximation Equation (17) in Equation (16), we obtain Δ , which, in conjunction with Green's theorem and the no-penetration boundary conditions on the solid walls, leads to the following Equation (19) is rewritten in the following matrix form , or , where , ∇ ∇ , are, respectively, the coefficients of matrix and . It must be noticed that is an eigenvalue of matrix . In this way the problem is solved by calculating the matrix C eigenvalues. The numerical solver takes as input a mesh, builds the matrix , and calculates its first eigenvalues and corresponding eigenvectors.

Validation of the Solver
This numerical solver has been tested first using meshes representing tanks filled with water but without cylinders (τ = 1) in the case of rectangular tank with vertical walls. In this case, the wavelength is twice the tank length for the first eigenmode. Consequently, for a given tank length , we have . Calculations were run with a 1.2-m-long tank and a width of 0.3 m. Three different levels of refinement were tested as shown in Figure 6, and the result concerning the computed field is presented in Figure 7. Figure 8 illustrates the convergence of the results concerning the wave number of the first eigenmode, and the corresponding numerical values are reported in Table 1.

Numerical Results against the Experiment
After the validation of the FEM solver for rectangular tank with vertical walls against the theoretical value of k for a tank without cylinders, the present method was tested against the experimental value for a dense array of cylinders presented in Molin et al. [26]; see Figure 9. In the above experiments only the eigenfrequency was reported, whereas the numerical FEM solver provides the wave number. In order to compare both results, an eigenfrequency is associated to the wave number using the dispersion relation of linearized water waves where ℎ 0.235 m is the water depth of the experimental tank (kh = 0.615); see Figures 1 and 9a. Tables 2 and 3 associate wave numbers from the solver and the eigenfrequencies due to Equation (10). The results were compared for different numbers of cylinders from 48 to 120. For the first comparison presented in Table 2 the meshes contained an average of 2748 nodes. Then a second comparison has been made with finer meshes with an average of 10,254 nodes, see Table 3.  Table 3. Numerical eigenfrequencies with an average of 10,254 nodes. The results were very encouraging. Indeed, on one hand numerical results were close to the experiment with an error systematically under 4.5%; see Table 4. On the other hand, increasing the level of refinement improved the accuracy of the results, see Figure 10. Unfortunately, this arrangement with cylinders very close to the tank wall was costly in computation time, and no better meshes were tested.

Independence of against the Number of Cylinders
A first set of calculations concerns the relation between the number of cylinders and the wave number for a given value of porosity. Three different meshes were tested with respectively 4, 36 and 100 cylinders, as shown in Figure 11. The quantities that remained the same in all variants were


The porosity τ;  Length and width;  Nx/Ny ratio. Porosity values of τ = 0.3, 0.5, 0.7 and 0.9 were tested and results are listed in Table 5. For all cases, the influence of the number of cylinders was found to be negligible. As a consequence, further calculations were undertaken with a number of cylinders as small as possible in order to reduce the computation time. Indeed, meshes with a small number of cylinders require a smaller number of nodes to coincide with all the details of the mesh.

Scaling Effect
A second set of calculations concerns the study of scaling effects. An initial pattern corresponding to the misaligned arrangement was repeated until reaching a tank length of L = 1.8 m; see Figure 12. The objective was to quantify the deviation of the non-dimensional wave number when the pattern was repeated. Obviously, this quantity did not remain the same and the results listed in Table  6 demonstrated that the deviation was not very important except for small porosities. Nevertheless, the deviation to the mean for small value of porosity remained a good approximation for general studies. Consequently, further computations are run with the smallest repeatable pattern.

FEM Convergence and Validation
In order to extend the previous results and study their property of repeatability, a series of numerical tests was performed for a tank size of L = 0.6 m and W = 0.3 m containing two cylinders and using both aligned and misaligned arrangements. One hundred meshes were generated with the mesh generator presented in Section 2.2 for porosity value ranging from τ = 0.3 up to τ = 0.95; see Figures 13 and 14. In Figure 15, the calculated results were plotted for two different levels of refinement. For the whole range of porosity value from 0.3 to 0.95, no bigger difference than 0.23% for the misaligned arrangement and 0.63% for the aligned arrangement was observed between the two levels of refinement. Thus, it was assumed that the average of 4186 nodes was enough to exploit the results.

Simplified Dispersion Equation in the Studied Porous Medium
From the above set of computations, Figure 16 was plotted for both aligned and misaligned arrangements. In the same figure the analytical expression provided by Molin et al. [26], Equation (16) that was considered among other formulations was also plotted for comparison. Figure 16. Results of the numerical campaign concerning G(τ) shown by dots. The expression G(τ) = τ 1/2 considered among other formulations by Molin et al. [26] is also shown by using the solid line.
Several polynomial expressions of G(τ) were proposed based on a least-square fit, as obtained using the constrained G(τ) = 1 for τ = 1. For the aligned case we obtained from the regression analysis for polynomial degree 2, 0.0508 0.5474 0.5033, 0.3 0.95, Furthermore, for > 0.5, the curve looks linear enough so that a linear regression was also proposed as follows: 0.4651 0.5349, 0.5 0.95.
The error between the expression of G(τ) and the numerically tested cases is presented in Table  7 below: The error between the expression of G(τ) and the numerical case is presented in the Table 8 below: which is the average of Equations (26) and (28) and was plotted in Figure 19. We got very good results from this regression, and the deviation did not become larger than 1.31%.  In concluding this part we could observe that both the present results agreed well with the dispersion relation proposed by Molin et al. [26] for porosities above 0.7. If for porosity down to 0.3, a better fit was obtained through the polynomial approximation of degree 3, and for porosity above 0.5, a linear fit appeared appropriate.

Estimation of in the Case of Random Arrangements
It was demonstrated in Section 3.3 that the dispersion equation could be very well modeled in the case of regular arrangements of cylinders. Nevertheless, in nature the arrangement is usually random. In this section, how relevant it is to model the studied porous medium based on random arrangements and extending the previous expressions was check. Ten random arrangements were tested with 16 cylinders per mesh for given porosity value of τ = 0.6, 0.7, 0.8 and 0.9 and were compared with the expression given by Equation (29). Results concerning G(τ) are shown in Figure  20 and presented a deviation in the random arrangement cases as compared to the regular arrangement; see also Table 9. It is observed that the model tended to overestimate the nondimensional wave number in most cases, still, however, the derived expression was useful in describing the macroscopic effect of porosity on the porous medium, especially for larger values of porosity.  Concerning the extension of present results to more general configurations it is seen that the number of cylinders was not important and the most significant parameters were the porosity and the cylinder arrangement. The latter were required for a direct application of these results without concern about the scale. Furthermore, if the tank or domain length/width ratio is not restricted, the error is found to be acceptable in most cases, and the deviation is less important for high porosity values and large number of cylinders. In the present work the discussion on the effect of the porosity, the specific surface, through the number of cylinders, and the arrangement of the cylinders on the local wavenumber within the porous media was carried out through sloshing conditions within a tank of given length, but valid whatever the water depth. Through the eigenvalues of the problem of sloshing, the calculated fundamental (or 1st harmonic) wave frequency had a wavelength λp at the scale of the porous media, given by λp = 2Lp. Calculation of the eigenvalues for harmonic oscillations of order n gives the frequency corresponding to λp = 2Lp/n. A complete information on the relation between the wavelength and the frequency on a given frequency range can be achieved by means of the present numerical method for various tank lengths, and this study is planned for future work.

Discussion and Conclusions
Different modifications are proposed concerning the dispersion relation through a dense array of cylinders considered as a porous medium within an enclosed domain. The derived expressions were established using different regular and irregular arrangements. As already demonstrated  [24], it was observed that for a given porosity and for a given arrangement, the wavenumber did not depend on the number of cylinders, or in other words, on the specific surface, defined as the solid-fluid contact surface per volume unit. According to the results presented, the modification of the dispersion relation was strongly dependent to the arrangement. These results were consistent with the numerical results based on the integral matching method, which showed through an interference process that the wavelength depended on the arrangement as shown when staggered or aligned cylinders were considered (see Figure 17 of Rey et al. [30]). In the case of regular arrangements, the wavenumber evolved linearly to the value provided by the standard linear water-wave dispersion relation for porosity larger than 0.7. Then, for low porosity (lower than 0.5) the wavenumber in the porous medium suddenly fell away from the water-wave value and this finding was repeated for all tested configurations. The present results need to be validated before being applied to general arrangements and configurations, and this task, including investigation of a very high value of porosity that requires increased computation cost, is planned for future work. Finally, as concerns the sloshing problem in a tank with vertical walls and structures, if flow separation effects are considered weak (in case of slow motions), the linear approach provides a useful approximate solution. In this case the present FEM is a very efficient computational tool providing useful information concerning the natural sloshing modes applicable to more general shape of tanks, with arbitrary distributions of internal structures with general cross section, which could be further exploited in CFD techniques generalizing the solution to include non-linear and viscous flow effects.