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

Joffrey Jamain, Julien Touboul, Vincent Rey, Kostas Belibassakis

► To cite this version:

Joffrey Jamain, Julien Touboul, Vincent Rey, Kostas Belibassakis. Porosity Effects on the Dispersion Relation of Water Waves through Dense Array of Vertical Cylinders. Journal of Marine Science and Engineering, 2020, 8 (12), pp.960. 10.3390/jmse8120960 . hal-03114584

HAL Id: hal-03114584

https://hal.science/hal-03114584

Submitted on 20 Jan 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Article

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

Joffrey Jamain<fnref n="1" />, Julien Touboul<fnref n="1" /> ID, Vincent Rey<fnref n="1" /> and Kostas Belibassakis<fnref n="2" />,*

Received: 19 October 2020; Accepted: 21 November 2020; Published: 24 November 2020

<figure><img src="image_3.png" /></figure>

check for
updates

**Abstract:** 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.

**Keywords:** water waves; porous medium; dispersion relation; FEM analysis

## 1. 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,

<footnote n="1">¹ Université Toulon, Aix Marseille Université, CNRS/INSU, IRD, MIO UM 110, Mediterranean Institute of Oceanography, 83130 La Garde, France; joffrey.jamain.20@seatech.fr (J.J.); julien.touboul@mio.osupytheas.fr (J.T.); vincent.rey@univ-tln.fr (V.R.)</footnote>

<footnote n="2">² School of Naval Architecture and Marine Engineering, National Technical University of Athens, Zografos, 15773 Athens, Greece</footnote>

<footnote>* Correspondence: kbel@fluid.mech.ntua.gr; Tel.: +30-210-7721138</footnote>
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.

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Sloshing tank used by Molin et al. [26].</figcaption></figure>

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 $\tau$ 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 $k_\tau$ 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.

## 2. Mathematical Model

### 2.1. 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 **u** is obtained as the gradient of a 3D velocity potential **Φ**:

$$ \mathbf{u} = \nabla_3 \Phi. \qquad (1) $$

Equation (1) in conjunction with the incompressibility of the water, leads to the Laplace equation:

$$ \Delta_3 \Phi = 0. \qquad (2) $$

We consider the condition of impermeability at the bottom of the tank ($z = -h$):

$$ \left. \frac{\partial \Phi}{\partial z} \right|_{z = -h} = 0, \qquad (3) $$

and the linearized free surface boundary condition at $z = 0$:

$$ \left. \frac{\partial^2 \Phi}{\partial t^2} \right|_{z=0} + g \left. \frac{\partial \Phi}{\partial z} \right|_{z=0} = 0. \qquad (4) $$

In the case of harmonic motion at given frequency $\omega$, exploiting the fact that the walls of cylinders and tank boundaries are vertical the evanescent modes are not excited and the wave field $\Phi$ is represented by:

$$ \Phi = f(z)\phi(x, y)e^{-i\omega t}, \qquad (5) $$

where $f$ is the vertical profile of the potential field and $\phi$ the corresponding 2D velocity potential on the x-y plane. With this form of potential, Equation (2) leads to:

$$ \frac{f''(z)}{f(z)} = -\frac{\Delta\phi(x, y)}{\phi(x, y)}, \qquad (6) $$

where a prime denotes differentiation with respect to the involved variable and $\Delta$ stands for the two-dimensional Laplace operator. This means that the quantity $-\frac{\Delta\phi}{\phi}$ is constant everywhere in the water domain. The latter is denoted by $k_{\tau}^2$, and thus, $k_{\tau}^2$ represents the wave number on the horizontal plane. Consequently, Equation (6) is written:

$$ f'' = k_{\tau}^2 f. \qquad (7) $$

Furthermore, from Equation (7) we obtain that the 2D potential $\phi$ satisfies the Helmholtz equation on the horizontal plane:

$$ \Delta\phi = -k_{\tau}^2\phi, \qquad (8) $$

valid in the water domain.

The solution of Equation (7) satisfying the boundary condition (3) is:

$$ f(z) = \frac{\cosh(k_{\tau}(z+h))}{\cosh(k_{\tau}h)}. \qquad (9) $$
which substituted in the potential $\Phi(x, y, z, t) = f(z)\phi(x, y)e^{-i\omega t}$ and used in the free surface condition (4) leads to:

$$ \omega^2 = gk_{\tau}\tanh(k_{\tau}h). \qquad (10) $$

In the sequel $k_{\tau}$ is considered constant everywhere in the domain, defined as

$$ -k_{\tau}^2 = \frac{\Delta\phi}{\phi}. \qquad (11) $$

It is noticed here that the associated wavelength

$$ \lambda_{\tau} = \frac{2\pi}{k_{\tau}}, \qquad (12) $$

does not correspond to the water wavelength $\lambda_1$, which for the first sloshing mode is twice the tank length $L$:

$$ \lambda_1 = \frac{2\pi}{k_1} = 2L, \qquad (13) $$

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 $k_{\tau}$ and $k_1$

$$ G(\tau) = \frac{k_{\tau}}{k_1}, \qquad (14) $$

in order to associate the water wavelength to its frequency in the porous medium

$$ \omega^2 = gk_1G(\tau)\tanh(k_1G(\tau)h). \qquad (15) $$

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 $\tau = 0.3$ up to 0.95 and a novel formulation of $G(\tau)$ 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(\tau) = \tau^{1/2}$.

## 2.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 Timokha ([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 $R_{in}$ 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 Timokha ([27], Equation (4.185)). However, the present results are found in good agreement with analytical predictions obtained by the variational Trefftz method by Timokha ([28], Figure 2).
<figure><img src="image_5.png" /><figcaption>**Figure 2.** The case of a circular 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). (a) Mesh, (b) first symmetric and (c) antisymmetric mode, as calculated by the present FEM.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The lowest eigenvalue corresponding to the symmetric mode in the case of a circular container of unit radius for various values of the radius *Rin* of the inner circular pole centered at *x* = 0.5 m considered by Timohka ([28], Section 4) as calculated by the present method.</figcaption></figure>

### 2.2.1. 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".

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

<figure><img src="image_7.png" /><figcaption>Figure 4. Different arrangements and mesh, before refinement (a) experimental; (b) misaligned; (c) aligned and (d) random arrangement.</figcaption></figure>

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

<figure><img src="image_8.png" /><figcaption>Figure 5. A 16 × 4 misaligned mesh repeated lengthwise and widthwise. This creates a 48 × 12 misaligned mesh.</figcaption></figure>

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

$$\iint_{\Omega} \Delta \phi W_i = -k^2 \iint_{\Omega} \phi W_i, \qquad (16)$$
where $(W_i)_{i=1,N}$ are the linear pyramid functions equal to 1 for node number $i$, and 0 elsewhere and $\Omega$ is the water domain on the horizontal plane. Denote the approximation of the field as

$$ \hat{\phi} = \sum_{j=1}^{N} \phi_j W_j, \qquad (17) $$

where the nodal values in the form of a global vector are

$$ X = \begin{pmatrix} \phi_1 \\ \vdots \\ \phi_N \end{pmatrix}. \qquad (18) $$

Using the FEM approximation Equation (17) in Equation (16), we obtain

$$ \sum_{j=1}^{N} \left[ \phi_j \iint_{\Omega} (\Delta W_j) W_i \right] = -k^2 \sum_{j=1}^{N} \left[ \phi_j \iint_{\Omega} W_j W_i \right], $$

which, in conjunction with Green's theorem and the no-penetration boundary conditions on the solid walls, leads to the following

$$ \sum_{j=1}^{N} \left[ \phi_j \iint_{\Omega} (-\nabla W_j \nabla W_i) \right] = -k^2 \sum_{j=1}^{N} \left[ \phi_j \iint_{\Omega} W_j W_i \right]. \qquad (19) $$

Equation (19) is rewritten in the following matrix form

$$ AX = -k^2BX, \qquad (20) $$

or

$$ B^{-1}AX = -k^2X, \qquad (21) $$

where

$$ A_{i,j} = - \iint_{\Omega} \nabla W_j \nabla W_i, \quad B_{i,j} = \iint_{\Omega} W_j W_i, \qquad (22) $$

are, respectively, the coefficients of matrix A and B. It must be noticed that $-k^2$ is an eigenvalue of matrix C = B$^{-1}$A. In this way the problem is solved by calculating the matrix C eigenvalues. The numerical solver takes as input a mesh, builds the matrix C = B$^{-1}$A, and calculates its first eigenvalues and corresponding eigenvectors.

## 2.4. Validation of the Solver

This numerical solver has been tested first using meshes representing tanks filled with water but without cylinders ($\tau = 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 L, we have $k_1 = \pi/L$. 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.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Empty tank and mesh (a) with 61 nodes; (b) with 213 nodes (refined once) and (c) with 793 nodes (refined twice).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** First sloshing eigenmode for tank without cylinders calculated by FEM mesh with 213 nodes. The wave field $\phi$ is presented using a colorplot.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** First sloshing eigenmode for tank without cylinders. Convergence of the FEM solver (log-log scale).</figcaption></figure>

<table><caption>**Table 1.** Wave number of an empty tank and error against the theory.</caption><thead><tr><td>Number of Nodes</td><td>Wave Number (m<sup>-1</sup>)</td><td>Error Against Theory</td></tr></thead><tbody><tr><td>61</td><td>2.6239</td><td>2.3%</td></tr><tr><td>213</td><td>2.6198</td><td>0.07%</td></tr><tr><td>793</td><td>2.6185</td><td>0.02%</td></tr></tbody></table>

## 2.5. 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.
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Sloshing tank with an arrangement of vertical cylinders: (a) Experimental set-up and (b) numerical tank and FEM model mesh (6913 nodes).</figcaption></figure>

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

$$ \omega^2 = gk_{\tau}\tanh(k_{\tau}h), \qquad (23) $$

where $h = 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.

<table><caption>**Table 2.** Numerical eigenfrequencies with an average of 2748 nodes.</caption><thead><tr><td>Number of Cylinders</td><td>Numerical Wave Number (m<sup>-1</sup>)</td><td>Associated Eigenfrequency (rad/s)</td></tr></thead><tbody><tr><td>48</td><td>2.52</td><td>3.62</td></tr><tr><td>60</td><td>2.49</td><td>3.59</td></tr><tr><td>84</td><td>2.44</td><td>3.52</td></tr><tr><td>120</td><td>2.38</td><td>3.44</td></tr></tbody></table>

<table><caption>**Table 3.** Numerical eigenfrequencies with an average of 10,254 nodes.</caption><thead><tr><td>Number of Cylinders</td><td>Numerical Wave Number (m<sup>-1</sup>)</td><td>Associated Eigenfrequency (rad/s)</td></tr></thead><tbody><tr><td>48</td><td>2.49</td><td>3.59</td></tr><tr><td>60</td><td>2.46</td><td>3.55</td></tr><tr><td>84</td><td>2.39</td><td>3.46</td></tr><tr><td>120</td><td>2.32</td><td>3.37</td></tr></tbody></table>

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.

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.

<table><caption>**Table 4.** Numerical eigenfrequencies against the experiment presented in Molin et al. [26].</caption><thead><tr><td>Number of Cylinders</td><td>Experimental Eigenfrequencies (rad/s)</td><td>Difference with the Numerical Eigenfrequencies (Average of 2748 Nodes)</td><td>Difference with the Numerical Eigenfrequencies (Average of 10,254 Nodes)</td></tr></thead><tbody><tr><td>48</td><td>3.55</td><td>2.10%</td><td>1.21%</td></tr><tr><td>60</td><td>3.50</td><td>2.55%</td><td>1.39%</td></tr><tr><td>84</td><td>3.41</td><td>3.17%</td><td>1.49%</td></tr><tr><td>120</td><td>3.30</td><td>4.31%</td><td>2.08%</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Comparison between the numerical model and the experiment.</figcaption></figure>

### 3. Estimation of the Porosity Effect on the Dispersion Relation

#### 3.1. Independence of G(τ) 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 $\tau$;

* Length and width;

* $Nx/Ny$ ratio.

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Porosity $\tau = 0.5$ with (a) 4 cylinders; (b) 36 cylinders and (c) 100 cylinders in the tank.</figcaption></figure>

Porosity values of $\tau = 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.

<table><caption>**Table 5.** Effect of the number of cylinders.</caption><thead><tr><th>Porosity τ</th><th>Wave Number with 4 Cylinders (m<sup>-1</sup>)</th><th>Wave Number with 36 Cylinders (m<sup>-1</sup>)</th><th>Wave Number with 100 Cylinders (m<sup>-1</sup>)</th><th>Mean (m<sup>-1</sup>)</th><th>Max Deviation to the Mean</th></tr></thead><tbody><tr><td>0.3</td><td>1.6968</td><td>1.7110</td><td>1.7122</td><td>1.7067</td><td>0.58%</td></tr><tr><td>0.5</td><td>2.0772</td><td>2.0996</td><td>2.1009</td><td>2.0926</td><td>0.74%</td></tr><tr><td>0.7</td><td>2.2946</td><td>2.3030</td><td>2.3037</td><td>2.3004</td><td>0.25%</td></tr><tr><td>0.9</td><td>2.5033</td><td>2.5052</td><td>2.5052</td><td>2.5046</td><td>0.05%</td></tr></tbody></table>

#### 3.2. 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.
<figure><img src="image_15.png" /><figcaption>**Figure 12.** Porosity value $\tau = 0.5$ with a tank of length of (a) $L = 0.6$ m and 2 cylinders; (b) $L = 1.2$ m and 4 cylinders and (c) $L = 1.8$ m and 6 cylinders.</figcaption></figure>

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.

<table><caption>**Table 6.** Effect of the number of cylinders.</caption><thead><tr><th>Porosity τ</th><th>k<sub>τ</sub>/k<sub>1</sub> with<br>2 Cylinders (m<sup>-1</sup>)</th><th>k<sub>τ</sub>/k<sub>1</sub> with<br>4 Cylinders (m<sup>-1</sup>)</th><th>k<sub>τ</sub>/k<sub>1</sub> with<br>6 Cylinders (m<sup>-1</sup>)</th><th>Mean (m<sup>-1</sup>)</th><th>Max Deviation<br>to the Mean</th></tr></thead><tbody><tr><td>0.3</td><td>0.6153</td><td>0.6421</td><td>0.6464</td><td>0.6346</td><td>3.1%</td></tr><tr><td>0.5</td><td>0.7631</td><td>0.7832</td><td>0.7877</td><td>0.7780</td><td>1.9%</td></tr><tr><td>0.7</td><td>0.8523</td><td>0.8683</td><td>0.8700</td><td>0.8635</td><td>1.3%</td></tr><tr><td>0.9</td><td>0.9495</td><td>0.9527</td><td>0.9532</td><td>0.9518</td><td>0.2%</td></tr></tbody></table>

### 3.3. Estimation of G(τ) in the Case of Regular Arrangements

#### 3.3.1. 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 $\tau = 0.3$ up to $\tau = 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.

<figure><img src="image_16.png" /><figcaption>**Figure 13.** Different misaligned arrangements and meshes used for the numerical campaign with (a) $\tau = 0.3$; (b) $\tau = 0.6$; and (c) $\tau = 0.95$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 14.** Different aligned arrangements and meshes used for the numerical campaign with (a) $\tau = 0.3$; (b) $\tau = 0.6$; and (c) $\tau = 0.95$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 15. Calculated non-dimensional wavenumber for different levels of mesh refinement, misaligned arrangement.</figcaption></figure>

### 3.3.2. 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><img src="image_19.png" /><figcaption>Figure 16. Results of the numerical campaign concerning $G(\tau)$ shown by dots. The expression $G(\tau) = \tau^{1/2}$ considered among other formulations by Molin et al. [26] is also shown by using the solid line.</figcaption></figure>

Several polynomial expressions of $G(\tau)$ were proposed based on a least-square fit, as obtained using the constrained $G(\tau) = 1$ for $\tau = 1$. For the aligned case we obtained from the regression analysis for polynomial degree 2,

$$G(\tau) = -0.0508\tau^2 + 0.5474\tau + 0.5033, \quad 0.3 < \tau < 0.95, \qquad (24)$$

and for polynomial degree 3,

$$G(\tau) = 0.9484\tau^3 - 2.0054\tau^2 + 1.7998\tau + 0.2572, \quad 0.3 < \tau < 0.95. \qquad (25)$$
Furthermore, for $\tau > 0.5$, the curve looks linear enough so that a linear regression was also proposed as follows:

$$G(\tau) = 0.4651\tau + 0.5349, \quad 0.5 < \tau < 0.95. \qquad (26)$$

The error between the expression of $G(\tau)$ and the numerically tested cases is presented in Table 7 below:

<table><caption>**Table 7.** Comparison between different polynomial models for the aligned arrangement.</caption><thead><tr><th>Deg. of the Polynomial</th><th>Max. Deviation to the Model</th><th>Max. Deviation to the Model (%)</th></tr></thead><tbody><tr><td>1 ($\tau$ > 0.5)</td><td>0.0075</td><td>0.97%</td></tr><tr><td>2</td><td>0.0357</td><td>3.57%</td></tr><tr><td>3</td><td>0.0150</td><td>2.39%</td></tr></tbody></table>

In the case of the misaligned arrangement, the following expressions were obtained, respectively:

$$G(\tau) = -0.0586\tau^2 + 0.5822\tau + 0.4764, \quad 0.3 < \tau < 0.95, \qquad (27)$$

$$G(\tau) = 0.8044\tau^3 - 1.7165\tau^2 + 1.6445\tau + 0.2676, \quad 0.3 < \tau < 0.95 \qquad (28)$$

The error between the expression of $G(\tau)$ and the numerical case is presented in the Table 8 below:

<table><caption>**Table 8.** Comparison between different polynomial models for the misaligned arrangement.</caption><thead><tr><th>Deg.</th><th>Max. Deviation to the Mean</th><th>Max. Deviation to the Mean (%)</th></tr></thead><tbody><tr><td>1 ($\tau$ > 0.5)</td><td>0.0056</td><td>0.73%</td></tr><tr><td>2</td><td>0.0314</td><td>5.10%</td></tr><tr><td>3</td><td>0.0138</td><td>2.24%</td></tr></tbody></table>

Finally, for high porosity values, the two curves were very close, see Figures 17 and 18. This is not surprising because $k_\tau/k_1 \to 1$ when $\tau \to 1$. So, it was suggested to use the following expression

$$G(\tau) = 0.4763\tau + 0.5237, \quad 0.5 < \tau < 0.95, \qquad (29)$$

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

<figure><img src="image_20.png" /><figcaption>**Figure 17.** Polynomial regressions for $G(\tau)$ in case of aligned arrangements: (a) polynomial degree 1; (b) polynomial degree 2 and (c) polynomial degree 3.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Polynomial regressions for G(τ) in the case of a misaligned arrangement: (a) polynomial deg.1; (b) polynomial deg.2 and (c) polynomial deg.3.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Linear regression for τ > 0.5 for both aligned and misaligned arrangements.</figcaption></figure>

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.

### 3.4. Estimation of G(τ) 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 non-dimensional 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.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Validity test of the model against random cylinder arrangements.</figcaption></figure>

<table><caption>**Table 9.** Deviation from the model for non-determinist cases.</caption><thead><tr><td>Porosity τ</td><td>Max. Deviation to the Model</td><td>Max. Deviation to the Model (%)</td></tr></thead><tbody><tr><td>0.6</td><td>0.0914</td><td>12.72%</td></tr><tr><td>0.7</td><td>0.1139</td><td>15.33%</td></tr><tr><td>0.8</td><td>0.0773</td><td>9.34%</td></tr><tr><td>0.9</td><td>0.0785</td><td>8.98%</td></tr></tbody></table>

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 = 2L_p$. Calculation of the eigenvalues for harmonic oscillations of order n gives the frequency corresponding to $λ_p = 2L_p/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.

## 4. 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 experimentally by Arnaud et al. [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.

**Author Contributions:** Conceptualization, K.B., J.T., V.R.; Formal analysis, J.J.; Investigation, J.J.; Project administration, K.B.; Resources, K.B.; Software, J.J.; Supervision, K.B.; Validation, J.J., V.R.; Visualization, J.J.; Writing—original draft, J.J., J.T., Writing—review and editing, K.B., J.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank Angeliki Karperaki from the Lab. of Ship and Marine Hydrodynamics of NTUA for useful discussion concerning the implementation of the FEM code.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Mei, C.C. *The Applied Dynamics of Ocean Surface Waves*; World Scientific: Singapore, 1989.

2. Molin, B. *Hydrodynamique des Structures Offshore*; Technip: Paris, France, 2002. (In French)

3. Goda, Y. *Wave Forces on a Vertical Circular Cylinder: Experiments and a Proposed Method of Wave Force Computation*; PARI Report 000-08; PARI: Yokosuka, Japan, 1964.

4. John, F. On the motion of floating bodies II. Simple harmonic motions. *Commun. Pure Appl. Math.* **1950**, *3*, 45–101. [CrossRef]

5. Yang, C.; Ertekin, R.C. Numerical Simulation of Nonlinear Wave Diffraction by a Vertical Cylinder. *J. Offshore Mech. Arct. Eng.* **1992**, *114*, 36–44. [CrossRef]

6. Young, R.W. Added mass and damping of a vertical cylinder in finite-depth waters. *Appl. Ocean Res.* **1981**, *3*, 119–133. [CrossRef]

7. Sabunku, T.; Calisal, S. Hydrodynamic coefficients for vertical radiation of circular cylinders in finite depth. *Ocean Eng.* **1981**, *8*, 25–63. [CrossRef]

8. McNatt, J.C. Cylindrical Linear Water Waves and Their Application to the Wave-Body Problem. Ph.D. Thesis, The University of Edinburgh, Edinburgh, UK, 2016.

9. Wang, C.; Tay, Z. Very Large Floating Structures: Applications, Research and Development. *Procedia Eng.* **2011**, *14*, 62–72. [CrossRef]

10. Falnes, J. *Ocean Waves and Oscillating Systems. Linear Interactions Including Wave Energy Extraction*; Cambridge University Press: New York, NY, USA, 2004.

11. Balitsky, P.; Quartier, N.; Fernández, G.V.; Stratigaki, V.; Troch, P. Analyzing the Near-Field Effects and the Power Production of an Array of Heaving Cylindrical WECs and OSWECs Using a Coupled Hydrodynamic-PTO Model. *Energies* **2018**, *11*, 3489. [CrossRef]

12. Arnaud, G. *Houle à la côte Propagation, Impact et Ouvrages Innovants*; Université de Toulon: Toulon, France, 2016. (In French)

13. Massel, S.R.; Furukawa, K.; Brinkman, R.M. Surface wave propagation in mangrove forests. *Fluid Dyn. Res.* **1999**, *24*, 219–249. [CrossRef]

14. Mei, C.C.; Chan, I.-C.; Liu, P.L.-F.; Huang, Z.; Zhang, W. Long waves through emergent coastal vegetation. *J. Fluid Mech.* **2011**, *687*, 461–491. [CrossRef]
15. Linton, C.M.; Evans, D.V. The interaction of waves with arrays of vertical circular cylinders. *J. Fluid Mech.* **1990**, *215*, 549–569. [CrossRef]

16. Maniar, H.D.; Newman, J.N. Wave diffraction by a long array of cylinders. *J. Fluid Mech.* **1997**, *339*, 309–330. [CrossRef]

17. McIver, P. Water-wave propagation through an infinite array of cylindrical structures. *J. Fluid Mech.* **2000**, *424*, 101–125. [CrossRef]

18. Kagemoto, H.; Murai, M.; Saito, M.; Molin, B.; Malenica, Š. Experimental and theoretical analysis of the wave decay along a long array of vertical cylinders. *J. Fluid Mech.* **2002**, *456*, 113–135. [CrossRef]

19. Belibassakis, K.; Arnaud, G.; Rey, V.; Touboul, J. Propagation and scattering of waves by dense arrays of impenetrable cylinders in a waveguide. *Wave Motion* **2018**, *80*, 1–19. [CrossRef]

20. Fu, Z.-J.; Xie, Z.-Y.; Ji, S.-Y.; Tsai, C.-C.; Li, A.-L. Meshless generalized finite difference method for water wave interactions with multiple-bottom-seated-cylinder-array structures. *Ocean Eng.* **2020**, *195*, 1067–1086. [CrossRef]

21. Sollitt, C.K.; Cross, R.H. Wave transmission through permeable breakwaters. In Proceedings of the 13th Conference Coastal Engineering, Vancouver, BC, Canada, 10–14 July 1972; ASCE: New York, NY, USA, 1972; p. 1827–1846.

22. Madsen, O.S. Wave transmission through porous structures. *J. Waterw. Harb. Coast. Eng. Div. ASCE* **1974**, *100*, 169–188.

23. Yu, X.; Chwang, A.T. Wave Motion through Porous Structures. *J. Eng. Mech.* **1994**, *120*, 989–1008. [CrossRef]

24. Arnaud, G.; Rey, V.; Touboul, J.; Sous, D.; Molin, B.; Gouaud, F. Wave propagation through dense vertical cylinder arrays: Interference process and specific surface effects on damping. *Appl. Ocean Res.* **2017**, *65*, 229–237. [CrossRef]

25. Arnaud, G.; Rey, V.; Touboul, J.; Sous, D. Wave Propagation Through Dense Vertical Cylinder Arrays: 3D Experimental Study. *Int. J. Ocean Coast. Eng.* **2019**, *2*, [CrossRef]

26. Molin, B.; Remy, F.; Arnaud, G.; Rey, V.; Touboul, J.; Sous, D. On the dispersion equation for linear waves traveling through or over dense arrays of vertical cylinders. *Appl. Ocean Res.* **2016**, *61*, 148–155. [CrossRef]

27. Faltinsen, O.M.; Timokha, A.N. *Sloshing*; Cambridge University Press: New York, NY, USA, 2009.

28. Timokha, A. Analytically approximate natural sloshing modes and frequencies for an upright circular container with poles. *J. Eng. Math.* **2016**, *101*, 47–54. [CrossRef]

29. Johnson, C. *Numerical Solution of Partial Differential Equations by the Finite Element Method*; Cambridge University Press: New York, NY, USA, 1988.

30. Rey, V.; Arnaud, G.; Touboul, J.; Belibassakis, K. Water wave scattering by dense or sparse arrays of surface-piercing bodies by integral matching method. *Appl. Ocean Res.* **2018**, *75*, 132–142. [CrossRef]

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_24.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).