2 Formulation of the problem

The basis of the problem is to consider a rigid vertical cylinder with circular cross-section $S_{c}$ of radius $r=a$ extended from the seabed and piercing the free surface. We start with a Cartesian coordinate system $(x,y,z)$ with $z=0$ at the mean free surface. In a first step, we solve the semianalytical problem where the water depth $h(x,y)$ is assumed to be constant for $r>a$ , and we calculate the mean drift forces on the cylinder in the case of the theory of linearized water waves, as described below. In a second step, we place a set of vertical rigid inclusions with arbitrary cross-section in the domain $a<r\leqslant b$ around the cylinder and we allow the seabed to vary in this domain. We then use the MSE described in the following combined with an FEM Galerkin’s method to calculate the first-order free-surface elevation and the mean drift forces on the cylinder.

In the theory of linearized water waves and in the time harmonic regime, the first-order velocity potential ${\it\Phi}^{(1)}(x,y,z,t)=\text{Re}\{{\it\phi}^{(1)}(x,y,z)\text{e}^{\text{i}{\it\omega}t}\}$ satisfies

where ${\it\omega}$ is the angular wave frequency, $g$ is the gravitational acceleration, $h(x,y)$ is the water depth and $\boldsymbol{{\rm\nabla}}=(\partial /\partial x,\partial /\partial y)^{\text{T}}$ , ${\it\phi}_{zz}^{(1)}=\partial ^{2}{\it\phi}^{(1)}/\partial z^{2}$ . For a wave propagating in the $x$ direction with amplitude $a_{0}$ and wavenumber $k$ linked to ${\it\omega}$ by the dispersion equation ${\it\omega}^{2}=gk\tanh (kh)$ , the expression of the potential is given by

or in the cylindrical coordinate system $(r,{\it\theta},z)$

where $\text{J}_{n}$ are the first-kind Bessel functions and ${\it\epsilon}_{n}=2-{\it\delta}_{n0}$ , where ${\it\delta}_{nm}$ is the Kronecker delta function. The corresponding first-order free-surface elevation is given by

With the vertical cylinder, the velocity potential ${\it\Phi}^{(1)}$ can be decomposed into an incident potential and a scattered potential, ${\it\Phi}^{(1)}={\it\Phi}_{I}^{(1)}+{\it\Phi}_{S}^{(1)}$ , with ${\it\Phi}_{I}^{(1)}$ given by (2.3). In addition to (2.1), ${\it\Phi}_{S}^{(1)}$ must satisfy the Neumann boundary condition on the cylinder surface,

with $\boldsymbol{n}_{\boldsymbol{s}}$ the outward unit normal on $S_{c}$ , and Sommerfeld’s radiation condition

It follows that the total potential can be expressed as (see MacCamy & Fuchs Reference MacCamy and Fuchs1954)

with $\text{H}_{n}$ the Hankel functions and $^{\prime }$ the derivative with respect to the argument. The forces on a fixed body are calculated by the integration of the hydrodynamic pressure on the body surface,

The term $\boldsymbol{F}(t)$ can be decomposed into the sum of the first-order term $\boldsymbol{F}^{(\mathbf{1})}(t)$ and the second-order term $\boldsymbol{F}^{(\mathbf{2})}(t)$ . The second-order term, which is proportional to the square of the incoming wave amplitude, contains a time-independent component (the drift mean force $F_{m}^{(2)}$ ) and a fluctuating component at frequency $2{\it\omega}$ . The term $F_{m}^{(2)}$ is given by

where $\text{d}l$ is the infinitesimal element along the water line at $z=0$ . In our problem, we want to study a refraction–diffraction problem of linear water waves with some rigid vertical object with a varying topography of the seabed. Thus, we choose to employ the mild-slope equation, which simplifies this linearized problem by assuming that the bed is locally flat, so that the velocity potential can be written as

which is valid only if $h(x,y)$ satisfies the mild-slope approximation,

A vertical averaging of the problem leads to the MSE derived by Berkhoff (Reference Berkhoff1972) and Smith & Sprinks (Reference Smith and Sprinks1975),

where the nonlinear terms have been neglected, $c_{1}=c_{p}c_{g}$ is the product of the phase velocity with the group velocity and $c_{2}=c_{g}/c_{p}$ is the ratio of the group and phase velocities.

Figure 1.

Schematic representation of the physical domains considered in the homogenization process. An elementary cell ${\it\zeta}Y_{i}^{\ast }$ containing an obstacle with Neumann data on its boundary is repeated periodically in the domain ${\it\Omega}_{{\it\zeta}}$ . When ${\it\zeta}\ll 1$ , one has a large number of small cells: the set of union of all elementary cells becomes dense in  ${\it\Omega}$ .

3 Homogenization theory applied to the MSE

We now use a geometric transformation in the MSE to design a cloak. Such an approach will introduce an anisotropic behaviour of the coefficient $c_{1}$ and a modification of the coefficient $c_{2}$ . The difficulty is then to represent this anisotropic behaviour of $c_{1}$ . For this, we propose to apply the theory of homogenization based on asymptotic methods (see Papanicolau, Bensoussan & Lions Reference Papanicolau, Bensoussan and Lions1978) to the MSE. The idea is then to replace the domain with anisotropic properties with a homogenized domain consisting of a periodic structure where the basic element is small with respect to the wavelength. In our case, the homogenized domain is called ${\it\Omega}_{{\it\zeta}}$ and the basic element ${\it\zeta}Y_{i}^{\ast }$ . The small positive parameter ${\it\zeta}$ acts as a scale factor on the unit elements $Y_{i}^{\ast }$ . A unit element $Y=[0,d_{1}]\times [0,d_{2}]$ , with a rigid vertical inclusion of arbitrary cross-section $S$ piercing the free surface, so that $Y^{\ast }=Y\setminus S$ , as seen in figure 1. The principal assumption of the problem is that the depth $h(x)$ ( $x=(x_{1},x_{2})$ ) is defined as a periodic function on $Y^{\ast }$ , with the same periodicity for $c_{p}(x)$ and $c_{g}(x)$ . Consequently, $h(x)$ , $c_{p}(x)$ and $c_{g}(x)$ on ${\it\Omega}_{{\it\zeta}}$ are now assumed to depend on a rapidly varying space variable $\tilde{x}=x/{\it\zeta}$ , so that $h(x)\rightarrow h(\tilde{x})$ , $c_{p}(x)\rightarrow c_{p}(\tilde{x})$ and $c_{g}(x)\rightarrow c_{g}(\tilde{x})$ . The periodicity of these functions is ${\it\zeta}d_{1}$ , ${\it\zeta}d_{2}$ with respect to the directions $\tilde{x_{1}}$ , $\tilde{x_{2}}$ . The condition on the boundary of the inclusion $\partial S_{i}$ is the no-flow condition. Altogether, the problem in ${\it\Omega}_{{\it\zeta}}$ in the mild-slope approximation can be written as

with $\partial S=\bigcup _{i}\partial S_{i}$ . The use of Einstein’s notation in $(\mathscr{P}_{{\it\zeta}})$ leads to

Defining ${\it\phi}$ and ${\it\sigma}_{i}({\it\phi})$ as asymptotic expansions with respect to ${\it\zeta}$ ,

it follows that

and

and similar expressions for higher-order terms of ${\it\sigma}_{i}$ . With consideration of the previous expressions, the problem $(\mathscr{P}_{{\it\zeta}})$ leads to

Therefore, identifying the terms that are factors of the same power of ${\it\zeta}$ , we find the following two equations:

Equation (3.7a ) is known to represent the behaviour of the studied physical quantity (velocity potential in the present case) at small scale in the periodic structure (annex problem A), and (3.7b ) gives the global behaviour of the physical quantity in the presence of the total structure (homogenized problem H). To calculate ${\it\phi}^{(0)}$ and ${\it\phi}^{(1)}$ appearing in (3.7) through the ${\it\sigma}^{i}$ , one requires the introduction of the mean operator $\langle \cdot \rangle$ acting on a function $f(\tilde{x})$ on $Y_{i}^{\ast }$ . By applying the mean operator to $(\text{H})$ ,

where the commutative property of $\langle \cdot \rangle$ with $\partial /\partial x_{i}$ has been used. The assumption of periodic functions on $Y_{i}^{\ast }$ and the no-flow condition on $\partial S_{i}$ leads to

where the divergence theorem was used. The homogenized problem $(\text{H})$ depends only on the space variable $x$ , which confirms that this equation represents the global effect of the structure. On the other hand, the annex problem $(\text{A})$ can be written as

where (3.5) was used. This expression of $(\text{A})$ can be interpreted as an equation for ${\it\phi}^{(1)}$ with the same periodicity as $h(\tilde{x})$ and parametrized by $x$ . The parameter $x$ appears in (3.10) in the derivative of ${\it\phi}^{(0)}$ only, so that with linearity ${\it\phi}^{(1)}$ can be expressed as

where the functions $w^{j}(\tilde{x})$ ( $j=1,\,2$ ) are solutions of (3.10) and correspond respectively to $\partial {\it\phi}^{(0)}/\partial x_{j}$ ( $j=1,\,2$ ) equal to $1$ . That is, $w^{j}(\tilde{x})$ are solutions of

and are periodic functions with respect to $\tilde{x}$ with period ( ${\it\zeta}d_{1}$ , ${\it\zeta}d_{2}$ ). Equation (3.12) depends on $\tilde{x}$ only and especially through $c_{1}(\tilde{x})$ . Therefore, $w^{j}$ can be calculated once and for all, regardless of ${\it\Omega}_{{\it\zeta}}$ . Given $w^{j}$ and applying the mean operator, (3.5) becomes

with

or

One can note that $\langle c_{1}\partial _{1}w^{2}\rangle =\langle c_{1}\partial _{2}w^{1}\rangle$ (see Guenneau, Zolla & Nicolet Reference Guenneau, Zolla and Nicolet2007). The problem (3.1) at zeroth order represents the so-called homogenized problem:

It should be noted that if the coefficients $c_{1}$ and $c_{2}$ in (3.1) are functions of both slow $(x)$ and fast $(\tilde{x})$ variables, then the homogenized system has homogenized coefficients depending on the slow variable $(x)$ (see Guenneau et al. Reference Guenneau, Zolla and Nicolet2007).

4 Cloaking with geometric transformation

We now follow the concept introduced by Pendry et al. (Reference Pendry, Schurig and Smith2006), and we apply it to the MSE for the design of a cylindrical cloak around a vertical cylinder. This leads to non-physical considerations over the bathymetry, and also with respect to the mild-slope approximation. Thus, as a second step, we discuss and numerically adapt the parameters resulting from the geometric transformation to our problem. Then, we use homogenization theory to represent the parameters by a set of rigid vertical inclusions with a specific layout. As noted before, we only consider the problem in its eigenbasis $(r,{\it\theta})$ , so we use a conformal map to obtain the final structure in the Cartesian basis $(x,y)$ . The selected geometric transformation is

with ${\it\alpha}=(b-a)/(b-r_{0})$ , ${\it\beta}=b\,(a-r_{0})/(b-r_{0})$ , $a$ the inner cylinder radius, $b$ the cloak outer radius and a small parameter $r_{0}\ll a$ introduced to avoid singularity of the transform at $r=a$ as in Kohn et al. (Reference Kohn, Shen, Vogelius and Weinstein2008). Such a transform in the MSE leads to

where $\unicode[STIX]{x1D645}$ is the Jacobian matrix of the transformation, $T=\unicode[STIX]{x1D645}^{-1}\unicode[STIX]{x1D645}^{-T}\det (\unicode[STIX]{x1D645})=\unicode[STIX]{x1D64D}({\it\theta})\unicode[STIX]{x1D63C}\unicode[STIX]{x1D64D}(-{\it\theta})$ , with $\unicode[STIX]{x1D63C}$ a diagonal matrix of the eigenvalues of the transformation in the basis $(r,{\it\theta})$ and $\unicode[STIX]{x1D64D}({\it\theta})$ the in-plane rotation matrix. The expressions of $\unicode[STIX]{x1D63C}$ and $\det (\unicode[STIX]{x1D645})$ are

In our study, we fixed $a=1~\text{m}$ , $b=3~\text{m}$ , $r_{0}=0.1~\text{m}$ , and the corresponding ${\it\lambda}_{i}$ are shown in figure 2(a). The outer radius $b$ can be as large as we want, but we choose $b=3~\text{m}$ to have a sufficient and realistic size of the structure. At this point, we need to modify the ${\it\lambda}_{i}$ in order to match the physical constraints.

Figure 2.

(a) Parameters ${\it\lambda}_{i}$ of the geometric transformation for $a=1~\text{m}$ , $b=3~\text{m}$ , $r_{0}=0.1~\text{m}$ . (b) Parameters ${\it\lambda}_{i}^{\prime }$ defined to adapt the seabed at $r=b$ and parameters $\tilde{{\it\lambda}_{i}}$ .

At $r=b$ , we have $c_{2}=c_{g}/c_{p}$ outside the cloak and $\tilde{c_{2}}={\it\lambda}_{3}c_{g}/c_{p}$ inside the cloak, with ${\it\lambda}_{3}>1$ . This value of ${\it\lambda}_{3}$ means that the seabed is discontinuous at $r=b$ . To avoid this discontinuity, we normalize ${\it\lambda}_{3}(r)$ by ${\it\lambda}_{3}(r=b)$ and accordingly, as ${\it\lambda}_{3}(r=b)={\it\lambda}_{2}(r=b)$ , we also normalize ${\it\lambda}_{2}(r)$ by ${\it\lambda}_{3}(r=b)$ . These new parameters are called respectively ${\it\lambda}_{3}^{\prime }$ and ${\it\lambda}_{2}^{\prime }$ and are shown in figure 2(b).

On the other hand, the coefficient ${\it\lambda}_{2}^{\prime }$ appears as a factor of $c_{1}=c_{p}c_{g}$ , and we know that after a certain threshold on $kh$ , the influence of variation of bathymetry becomes negligible. Therefore, for convenience the value of ${\it\lambda}_{2}^{\prime }(r=a)$ is fixed. Then the curve $\tilde{{\it\lambda}}_{2}$ presented in figure 2(b) is obtained by matching ${\it\lambda}_{2}^{\prime }$ with the value of ${\it\lambda}_{2}^{\prime }(r=a)$ .

The last consideration concerns the value ${\it\lambda}_{1}$ at $(r=a)$ . The term ${\it\lambda}_{1}$ represents the anisotropy in the radial direction. In the sense of homogenization, a small value of ${\it\lambda}_{1}$ means a small thickness of the inclusion in the radial direction, which can be non-realistic. We thus need to increase the value ${\it\lambda}_{1}(r=a)$ in order to have a sufficient thickness of the inclusion with regard to wave loads. The curve $\tilde{{\it\lambda}}_{1}$ is defined in the same manner as $\tilde{{\it\lambda}_{2}}$ (see figure 2 b). Accordingly, the term $\tilde{c}_{1}$ is now defined as $\tilde{c}_{1}=\text{diag}(\tilde{{\it\lambda}}_{1},\tilde{{\it\lambda}}_{2})\,c_{1}$ .

Figure 3.

Results of homogenization and conformal mapping. (a) Discrete bathymetry (dashed grey curve) resulting from homogenization calculus and its continuous approximation (solid black line) for $a\leqslant r\leqslant b$ . The top view of (a) represents the repartition of inclusions before conformal mapping. (b) Result of the conformal mapping for the domain in (a).

Meanwhile, $\tilde{{\it\lambda}}_{3}$ is a result of the homogenization process. The adaptation of the parameters given by the geometric transform has been carried out for $k_{0}h_{0}=0.55$ , meaning that we work in intermediate/shallow water. We can now perform homogenization in coordinates $(r,{\it\theta})$ , and use the conformal map defined in (4.4) to obtain the final circular structure in $(x,y)$ coordinates. The conformal map that maps a rectangular domain ${\it\Omega}_{1}:\;[a,b]\times [-{\rm\pi}/{\it\gamma},{\rm\pi}/{\it\gamma}]$ on an annular domain ${\it\Omega}_{2}:\;[a,b]\times [0,2{\rm\pi}]$ is given by

where $z$ is the complex coordinate in the domain ${\it\Omega}_{1}$ . The transformation (4.4) allows us to solve the annex problem (3.12) on a rectangular cell and transpose the result in the domain of interest $a\leqslant r\leqslant b$ around the cylinder, keeping all the metric properties (see figure 3 b). Homogenization calculus with respect to the conformal map defined above requires the discretization of $\tilde{{\it\lambda}}_{1}(r)$ and $\tilde{{\it\lambda}}_{2}(r)$ on a number $N_{L}$ of layers, where the thickness of each layer depends on the number of angular sectors we want in them, as shown in figure 3(a). It is important to note that we have the larger layers where the variations of the different parameters are strong. This is due to the fact that the conformal map conserves the surface of the inclusions, which implies that the thickness becomes very small. In addition, we solve (3.12) for all of the layers on a unit cell $Y:\,[0.2\times 1]$ with an inclusion of rectangular cross-section and we verify the equality between (3.15) and $\tilde{c}_{1}$ . For each layer of a given water depth, the homogenization gives the dimensions of the inclusions (see figure 3 a). Then, the cells $Y$ are scaled with respect to the thickness of each layer and mapped onto the annular domain around the cylinder. These results are summarized in figure 3 and table 1 (structure 3), where we give the radial thickness and angular sector of the inclusions in each layer. As we find a discrete bathymetry with this process, we interpolate it with respect to the mild-slope approximation as shown in figure 3(a), and the corresponding values for $\tilde{c_{2}}$ are represented in figure 2(b) through $\tilde{{\it\lambda}}_{3}(r)$ .

Table 1.

Parameters for the different structures studied. The column entitled layer represents the numbering of the layers starting from $r=a$ , the $N_{i}$ columns are the number of inclusions per layer, the ${\rm\Delta}{\it\theta}$ columns are the angular sector in radians of the inclusions in each layer and the ${\rm\Delta}R$ columns are the radial thickness in centimetres of the inclusions in each layer.

5 Numerical results

The result of the process described above to construct the structure is not unique. For the fixed bathymetry presented in figure 3 (solid black curve) we find other structures by changing the number of layers and the number of inclusions in each layer. We have tested several structures and we choose to present four structures, for clarity, where the parameters are given in table 1. Because the water depth is varying more rapidly towards the interior of the cloak, the mild-slope approximation is only valid in the following for $kh>0.2$ . We are interested in two arguments to characterize the hydrodynamic behaviour of the devices. The first one is a quantity called a ‘cloaking factor’, which is not defined as in Porter & Newman (Reference Porter and Newman2014), but as follows. We integrate the norm of the scattered potential from $r=b+3{\it\lambda}$ to $r=b+5{\it\lambda}$ over ${\it\theta}=[0;2{\rm\pi}]$ with ${\it\lambda}$ the wavelength. We call $\text{SF}_{c}$ the case with a constant flat seabed and $\text{SF}_{s}$ the case with the cylinder surrounded by the structure, and so the cloaking factor appears as

With this definition, when $\text{CF}<1$ , the cylinder surrounded by the structure scatters less energy than with a flat bed, and perfect cloaking is reached for $\text{CF}=0$ . An important fact in this study is the successive approximations made in § 4 in conjunction with the parameters ${\it\lambda}_{i}^{\prime }$ and $\tilde{{\it\lambda}}_{i}$ . Each of these defines an anisotropic fluid that departs from perfect cloaking. Therefore, by performing calculations of the cloaking factor for all of the anisotropic fluids that we defined, we can quantify how far we are from perfect cloaking before the homogenization process. This gives an idea of the error introduced in the ideal model of cloaking and of the error due to homogenization. Results of these calculations are given in figure 4(a). The first approximation, which is the use of the parameter $r_{0}$ , leads to a cloaking factor close to zero for all of the $kh$ range. The second and third approximations, related respectively to the ${\it\lambda}_{i}^{\prime }$ and $\tilde{{\it\lambda}}_{i}$ , make the model deviate more and more from perfect cloaking but in an acceptable way for $kh\sim 0.5$ . Finally, the homogenization process performed from the $\tilde{{\it\lambda}}_{i}$ provides us with a good approximation of the anisotropic fluid, as shown by the similarity of the curve trends.

Figure 4.

(a) The cloaking factor $\text{CF}$ as defined in (5.1) against $kh$ for the four structures. The curves denoted as $r_{0}$ , ${\it\lambda}^{\prime }$ and $\tilde{{\it\lambda}}$ correspond to calculations of the cloaking factor for the anisotropic fluids defined at each step of the successive approximations detailed in § 4. (b) The mean drift force $F_{m}^{(2)}$ (see (2.9)) against $kh$ for the four structures (dotted lines) compared with the case of a cylinder alone on a constant flat bed (solid line).

The second argument is the mean drift force given by (2.9). The results are normalized by ${\it\rho}ga_{0}^{2}a$ , and we compare the case with a flat bed and the case with the structure. Figure 4(a) gives the cloaking factor with respect to $kh$ for the four structures tested. Except for structure number 4, $\text{CF}$ has a quasisimilar form to the others. We can see a reduction of the scattered energy of greater than 10 % for $kh\in [0.2;0.8]$ ; this interval can be extended to $kh$ close to 1 for some structures, such as number 3. On the other hand, the scattered energy is divided by a factor greater than 2 for $kh\in [0.25;0.55]$ , which corresponds to a structure placed in the coastal zone. We obtain the best cloaking effect in the previous interval with a reduction of the scattered energy equal to 80 % for $kh=0.46$ (structure 3). Perfect cloaking has not been achieved in this study, but we have reduced the scattered energy for a range of $kh$ in shallow water and intermediate water. Because of the drastic modifications of the parameters obtained with geometric transform, it seems clear that we cannot achieve perfect cloaking but only approach it.

The introduction of a region around the cylinder with varying bathymetry and containing a set of rigid vertical inclusions modifies not only the scattered field but also the forces acting on the cylinder. Figure 4(b) represents the mean drift force $F_{m}^{(2)}$ with respect to $kh$ and the comparison with the case of a flat bed. We observe a reduction of $F_{m}^{(2)}$ for $kh\in [0.2;0.7]$ , which is included in the range of $kh$ where we reduced the scattered field. The best range of $kh$ appears to be $kh\in [0.2;0.6]$ , where we divide by two the mean drift force acting on the cylinder (except for structure 4). Searching for a device for which we reduce $F_{m}^{(2)}$ for $kh>1$ with this method leads to structure 4, but the cloaking factor is not so good. With the method we chose, i.e. the use of a simple geometric transformation to construct a cloaking device, the adaptation of the parameters and the use of homogenization theory, we found a device that worked quite well for $kh\in [0.2;0.7]$ for linear water waves. Obviously, a device with some rigid vertical objects piercing the free surface is probably strongly dissipative, and a more complete study would be needed. On the other hand, a numerical optimization of the aspect ratio of the inclusions and of the layers should lead to an improvement of the results.