1 Introduction
Ten years ago, a word emerged from the optical community and has grown so much in importance that it is now employed by all of the metamaterials community to denominate the research on devices that render an object invisible for far-field observers, namely ‘cloaking’. Pendry, Schurig & Smith (Reference Pendry, Schurig and Smith2006) proposed to use a change of coordinates in Maxwell’s equations, keeping them invariant, to define a coating surrounding an object to render it invisible, and Leonhardt (Reference Leonhardt2006), meanwhile, proposed to use conformal mappings in Helmholtz’s equation. The former approach results in a material with anisotropic physical parameters (permittivity and permeability in optics), where waves propagating from free space are bent around an invisibility area and appear undisturbed to any far-field observer. Unfortunately, this method introduces severe limitations in the design of cloaks, due to the singular behaviour of the material parameters at the inner boundary of the cloak. Since the introduction of these concepts, some work has been carried out to complete and extend (see Kadic et al. Reference Kadic, Buckmann, Schittny and Wegener2013) the first ideas to other wave systems such as in acoustics (see Norris Reference Norris2008) and elastodynamics (see Milton, Briane & Willis Reference Milton, Briane and Willis2006; Farhat et al. Reference Farhat, Guenneau, Enoch and Movchan2009).
Regarding water waves, the first study was performed by Farhat et al. (Reference Farhat, Enoch, Guenneau and Movchan2008), where a geometric transform was used to design a circular cloak with homogenization theory to represent the modifications of the wave equation by an array of vertical rigid inclusions piercing the free surface. Since then, some studies have been carried out in order to control water wave propagation. Alam (Reference Alam2012) considered stratified seas to cloak floating objects; Berraquero et al. (Reference Berraquero, Maurel, Petitjeans and Pagneux2013) realized a water wave shifter; Newman (Reference Newman2014) used additional bodies with optimized geometries surrounding a cylinder to give numerical evidence for perfect cloaking; McIver (Reference McIver2014) proved that perfect cloaking cannot be achieved with fixed bodies, floating bodies or variations in the seabed topography; Porter & Newman (Reference Porter and Newman2014) proposed a numerical optimization process over bathymetry around a circular cylinder in order to achieve near-perfect cloaking for linear waves; Dupont et al. (Reference Dupont, Kimmoun, Molin, Guenneau and Enoch2015) studied an invisibility carpet in a wavetank for linear waves; recently, Zareei & Alam (Reference Zareei and Alam2015) used nonlinear transformation to achieve cloaking in shallow-water waves.
Here, we propose to apply a linear geometric transform in the mild-slope equation (MSE) to cloak a rigid vertical cylinder. We choose the MSE to allow for a slow variation of the seabed, combined with the interaction of waves with rigid vertical objects. The use of a geometric transform induces physical limitations at the inner and outer boundaries of the cloaked region, so we propose first to adapt the parameters of our problem, and then to use homogenization theory to achieve these parameters. The final cloaking device is constructed using a conformal map. We perform numerical investigation using the finite element method (FEM) to show that we can approach a cloaking effect with non-ideal parameters and a reduction of the mean drift force (see Porter & Newman Reference Porter and Newman2014).
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
$$\begin{eqnarray}\displaystyle & \displaystyle {\rm\nabla}^{2}{\it\phi}^{(1)}+{\it\phi}_{zz}^{(1)}=0\quad \text{on }-h(x,y)<z<0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle g{\it\phi}_{z}^{(1)}-{\it\omega}^{2}{\it\phi}^{(1)}=0\quad \text{on }z=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\phi}_{z}^{(1)}+\boldsymbol{{\rm\nabla}}{\it\phi}^{(1)}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}h=0\quad \text{on }z=-h(x,y), & \displaystyle\end{eqnarray}$$
${\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 
$$\begin{eqnarray}{\it\phi}_{I}^{(1)}(x,y,z)=-\text{i}{\displaystyle \frac{a_{0}g}{{\it\omega}}}{\displaystyle \frac{\cosh (k(z+h))}{\cosh (kh)}}\text{e}^{\text{i}kx},\end{eqnarray}$$
or in the cylindrical coordinate system
$(r,{\it\theta},z)$
$$\begin{eqnarray}{\it\phi}_{I}^{(1)}(r,{\it\theta},z)=-\text{i}{\displaystyle \frac{a_{o}g}{{\it\omega}}}{\displaystyle \frac{\cosh (k(z+h))}{\cosh (kh)}}\displaystyle \mathop{\sum }_{n=0}^{\infty }{\it\epsilon}_{n}\text{i}^{n}\text{J}_{n}(kr)\cos (n{\it\theta}),\end{eqnarray}$$
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
$$\begin{eqnarray}\boldsymbol{{\it\eta}}_{I}^{(1)}=\text{Re}\{{\it\eta}_{I}^{(1)}\text{e}^{\text{i}{\it\omega}t}\},\quad \text{with }{\it\eta}_{I}^{(1)}=-{\displaystyle \frac{1}{g}}{\it\Phi}_{It}^{(1)}|_{z=0}.\end{eqnarray}$$
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,
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}{\it\Phi}_{S}^{(1)}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{ s}}=-\boldsymbol{{\rm\nabla}}{\it\Phi}_{I}^{(1)}\boldsymbol{\cdot }\boldsymbol{n}_{\boldsymbol{ s}}\quad \text{on }S_{c},\end{eqnarray}$$
with
$\boldsymbol{n}_{\boldsymbol{s}}$
the outward unit normal on
$S_{c}$
, and Sommerfeld’s radiation condition
$$\begin{eqnarray}\lim _{r\rightarrow \infty }\sqrt{r}\left({\displaystyle \frac{\partial {\it\Phi}_{S}^{(1)}}{\partial r}}+\text{i}k{\it\Phi}_{S}^{(1)}\right)=0.\end{eqnarray}$$
It follows that the total potential can be expressed as (see MacCamy & Fuchs Reference MacCamy and Fuchs1954)
$$\begin{eqnarray}{\it\phi}^{(1)}(r,{\it\theta},z)=-\text{i}{\displaystyle \frac{a_{o}g}{{\it\omega}}}{\displaystyle \frac{\cosh (k(z+h))}{\cosh (kh)}}\displaystyle \mathop{\sum }_{n=0}^{\infty }{\it\epsilon}_{n}\text{i}^{n}\left(\text{J}_{n}(kr)-{\displaystyle \frac{\text{J}_{n}^{\prime }(ka)}{\text{H}_{n}^{\prime }(ka)}}\text{H}_{n}(kr)\right)\cos (n{\it\theta}),\end{eqnarray}$$
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,
$$\begin{eqnarray}\boldsymbol{F}(t)=\int _{S_{c}}{\it\rho}\left({\it\Phi}_{t}+gz+{\displaystyle \frac{1}{2}}(\boldsymbol{{\rm\nabla}}{\it\Phi})^{2}\right)\boldsymbol{n}_{\boldsymbol{ s}}\,\text{d}S.\end{eqnarray}$$
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
$$\begin{eqnarray}F_{m}^{(2)}={\displaystyle \frac{1}{T}}\int _{0}^{\text{T}}\left[-{\displaystyle \frac{1}{2}}{\it\rho}g\int _{z=0}{{\it\eta}^{(1)}}^{2}\boldsymbol{n}_{\boldsymbol{ s}}\,\text{d}l+{\displaystyle \frac{1}{2}}{\it\rho}\int _{S_{c}}\boldsymbol{{\rm\nabla}}{\it\Phi}^{(1)2}\boldsymbol{n}_{\boldsymbol{ s}}\,\text{d}S\right]\,\text{d}t,\end{eqnarray}$$
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
$$\begin{eqnarray}{\it\phi}(x,y,z)\approx w(h,z){\it\phi}(x,y)\quad \text{with }w(h,z)={\displaystyle \frac{\cosh (k(z+h))}{\cosh (kh)}}\text{ and }k=k(h),\end{eqnarray}$$
which is valid only if
$h(x,y)$
satisfies the mild-slope approximation,
$$\begin{eqnarray}|\boldsymbol{{\rm\nabla}}h|/|kh|<O({\it\epsilon})\quad ({\it\epsilon}\ll 1).\end{eqnarray}$$
A vertical averaging of the problem leads to the MSE derived by Berkhoff (Reference Berkhoff1972) and Smith & Sprinks (Reference Smith and Sprinks1975),
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(c_{1}\boldsymbol{{\rm\nabla}}){\it\phi}+{\it\omega}^{2}c_{2}{\it\phi}=0,\end{eqnarray}$$
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
$$\begin{eqnarray}(\mathscr{P}_{{\it\zeta}}):\left\{\begin{array}{@{}ll@{}}-\displaystyle \mathop{\sum }_{i,j=1}^{2}{\displaystyle \frac{\partial }{\partial x_{i}}}\left(c_{1}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}(x)}{\partial x_{j}}}\right)=c_{2}(\tilde{x}){\it\omega}^{2}{\it\phi}(x)\quad & \text{in }{\it\Omega}_{{\it\zeta}},\\ c_{1}(\tilde{x})\displaystyle \mathop{\sum }_{i=1}^{2}{\displaystyle \frac{\partial {\it\phi}}{\partial x_{i}}}n_{i}=0\quad & \text{on }{\it\zeta}\partial S,\end{array}\right.\end{eqnarray}$$
with
$\partial S=\bigcup _{i}\partial S_{i}$
. The use of Einstein’s notation in
$(\mathscr{P}_{{\it\zeta}})$
leads to
$$\begin{eqnarray}-{\displaystyle \frac{\partial }{\partial x_{i}}}{\it\sigma}_{i}({\it\phi}(x))=c_{2}(\tilde{x}){\it\omega}^{2}{\it\phi}(x),\quad \text{with }{\it\sigma}_{i}({\it\phi}(x))=c_{1}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}(x)}{\partial x_{j}}}.\end{eqnarray}$$
Defining
${\it\phi}$
and
${\it\sigma}_{i}({\it\phi})$
as asymptotic expansions with respect to
${\it\zeta}$
,
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\it\phi}={\it\phi}^{(0)}(x)+{\it\zeta}{\it\phi}^{(1)}(x,\tilde{x})+{\it\zeta}^{2}{\it\phi}^{(2)}(x,\tilde{x})+O({\it\zeta}^{3}),\\ {\it\sigma}_{i}={\it\sigma}_{i}^{(0)}+{\it\zeta}{\it\sigma}_{i}^{(1)}+{\it\zeta}^{2}{\it\sigma}_{i}^{(2)}+O({\it\zeta}^{3}),\end{array}\right\}\end{eqnarray}$$
it follows that
$$\begin{eqnarray}{\displaystyle \frac{\partial {\it\phi}}{\partial x_{i}}}=\left({\displaystyle \frac{\partial {\it\phi}^{(0)}}{\partial x_{i}}}+{\displaystyle \frac{\partial {\it\phi}^{(1)}}{\partial \tilde{x}_{i}}}\right)+{\it\zeta}\left({\displaystyle \frac{\partial {\it\phi}^{(1)}}{\partial x_{i}}}+{\displaystyle \frac{\partial {\it\phi}^{(2)}}{\partial \tilde{x}_{i}}}\right)+O({\it\zeta}^{2})\end{eqnarray}$$
and
$$\begin{eqnarray}{\it\sigma}_{i}^{(j)}({\it\phi})=c_{1}(\tilde{x})\left({\displaystyle \frac{\partial {\it\phi}^{(j)}}{\partial x_{i}}}+{\displaystyle \frac{\partial {\it\phi}^{(j+1)}}{\partial \tilde{x}_{i}}}\right),\quad j=(0,1),\end{eqnarray}$$
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
$$\begin{eqnarray}-\!\left({\displaystyle \frac{1}{{\it\zeta}}}{\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}+{\displaystyle \frac{\partial }{\partial x_{i}}}\right)({\it\sigma}_{i}^{(0)}+{\it\zeta}{\it\sigma}_{i}^{(1)}+O({\it\zeta}^{2}))=c_{2}(\tilde{x}){\it\omega}^{2}({\it\phi}^{(0)}(x)+O({\it\zeta})).\end{eqnarray}$$
Therefore, identifying the terms that are factors of the same power of
${\it\zeta}$
, we find the following two equations:
$$\begin{eqnarray}\displaystyle & \displaystyle (\text{A}):\quad -{\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}{\it\sigma}_{i}^{(0)}=0\quad (\text{for }{\it\zeta}^{-1}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (\text{H}):\quad -{\displaystyle \frac{\partial }{\partial x_{i}}}{\it\sigma}_{i}^{(0)}-{\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}{\it\sigma}_{i}^{(1)}=c_{2}(\tilde{x}){\it\omega}^{2}{\it\phi}^{(0)}\quad (\text{for }{\it\zeta}). & \displaystyle\end{eqnarray}$$
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})$
,
$$\begin{eqnarray}(\text{H}):\quad -{\displaystyle \frac{\partial }{\partial x_{i}}}\langle {\it\sigma}_{i}^{(0)}\rangle -\left\langle {\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}{\it\sigma}_{i}^{(1)}\right\rangle ={\it\omega}^{2}\langle c_{2}(\tilde{x})\rangle {\it\phi}^{(0)},\end{eqnarray}$$
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
$$\begin{eqnarray}(\text{H}):\quad -{\displaystyle \frac{\partial }{\partial x_{i}}}\langle {\it\sigma}_{i}^{(0)}\rangle ={\it\omega}^{2}\langle c_{2}(\tilde{x})\rangle {\it\phi}^{(0)},\end{eqnarray}$$
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
$$\begin{eqnarray}(\text{A}):\quad -{\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}\left(c_{1}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}^{(1)}}{\partial \tilde{x}_{i}}}\right)=\left({\displaystyle \frac{\partial {\it\phi}^{(0)}}{\partial x_{i}}}\right)\left({\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}c_{1}(\tilde{x})\right),\end{eqnarray}$$
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
$$\begin{eqnarray}{\it\phi}^{(1)}(x,\tilde{x})={\displaystyle \frac{\partial {\it\phi}^{(0)}}{\partial x_{j}}}w^{j}(\tilde{x}),\end{eqnarray}$$
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
$$\begin{eqnarray}(\text{A}):\quad -{\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}\left(c_{1}(\tilde{x}){\displaystyle \frac{\partial w^{j}}{\partial \tilde{x}_{i}}}\right)={\it\delta}_{ij}\left({\displaystyle \frac{\partial }{\partial \tilde{x}_{i}}}c_{1}(\tilde{x})\right)\end{eqnarray}$$
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
$$\begin{eqnarray}\langle {\it\sigma}_{i}^{(0)}\rangle (x)=[c_{1}]_{ij}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}^{0}(x)}{\partial x_{j}}},\end{eqnarray}$$
with
$$\begin{eqnarray}[c_{1}]_{ij}(\tilde{x})={\displaystyle \frac{1}{\text{area}(Y^{\ast })}}\iint _{Y^{\ast }}c_{1}(\tilde{x})\left(1+{\displaystyle \frac{\partial w^{j}(\tilde{x})}{\partial \tilde{x}_{i}}}\right)\,\text{d}\tilde{x}\end{eqnarray}$$
or
$$\begin{eqnarray}[c_{1}](\tilde{x})=\left(\begin{array}{@{}cc@{}}\langle c_{1}\rangle -\langle c_{1}\partial _{1}w^{1}\rangle & -\langle c_{1}\partial _{1}w^{2}\rangle \\ -\langle c_{1}\partial _{2}w^{1}\rangle & \langle c_{1}\rangle -\langle c_{1}\partial _{2}w^{2}\rangle \end{array}\right).\end{eqnarray}$$
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:
$$\begin{eqnarray}(\mathscr{P}_{0}):\left\{\begin{array}{@{}ll@{}}-\displaystyle \mathop{\sum }_{i,j=1}^{2}{\displaystyle \frac{\partial }{\partial x_{i}}}\left([c_{1}]_{ij}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}^{(0)}(x)}{\partial x_{j}}}\right)={\it\omega}^{2}\langle c_{2}(\tilde{x})\rangle {\it\phi}^{0}(x)\quad & \text{in }{\it\Omega}_{{\it\zeta}},\\ \displaystyle \mathop{\sum }_{i,j=1}^{2}[c_{1}]_{ij}(\tilde{x}){\displaystyle \frac{\partial {\it\phi}^{(0)}}{\partial x_{j}}}n_{i}=0\quad & \text{in }{\it\zeta}\partial S.\end{array}\right.\end{eqnarray}$$
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
$$\begin{eqnarray}f:\left\{\begin{array}{@{}ll@{}}r^{\prime }={\it\alpha}\,r+{\it\beta}\quad & \text{on }0<r<b,\\ {\it\theta}^{\prime }={\it\theta}\quad & \text{on }0<{\it\theta}<2{\rm\pi},\end{array}\right.\end{eqnarray}$$
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
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(c_{1}T\boldsymbol{{\rm\nabla}}{\it\phi})+c_{2}\det (J){\it\omega}^{2}{\it\phi}=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(\tilde{c}_{1}\boldsymbol{{\rm\nabla}}{\it\phi})+\tilde{c_{2}}{\it\omega}^{2}{\it\phi}=0,\end{eqnarray}$$
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
$$\begin{eqnarray}\unicode[STIX]{x1D63C}=\text{diag}\,({\it\lambda}_{1},{\it\lambda}_{2})=\text{diag}\left({\displaystyle \frac{r-{\it\beta}}{r}},{\displaystyle \frac{r}{r-{\it\beta}}}\right),\quad \det (\unicode[STIX]{x1D645})={\it\lambda}_{3}={\displaystyle \frac{r-{\it\beta}}{{\it\alpha}^{2}r}}.\end{eqnarray}$$
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
$$\begin{eqnarray}w={\it\xi}\,\text{e}^{{\it\gamma}z},\quad \text{with }{\it\gamma}={\displaystyle \frac{\ln (b/a)}{b-a}},{\it\xi}={\displaystyle \frac{b+a}{\text{e}^{{\it\gamma}b}+\text{e}^{{\it\gamma}a}}},\end{eqnarray}$$
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
$$\begin{eqnarray}\text{CF}={\displaystyle \frac{\text{SF}_{s}}{\text{SF}_{c}}}.\end{eqnarray}$$
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.
6 Conclusion
We have conclusively shown using numerical simulations that it is possible to approximate a cloaking effect using a simplified set of parameters based on geometric transforms. In a first step, we considered a linear geometric transform applied to the MSE, and we modified this ideal result to respect the constraints at the limits of the cloak, i.e. we made the seabed continuous at the outer boundary of the cloak and we fixed a finite water depth at the inner boundary of the cloak. In a second step, we used homogenization theory combined with a conformal map to construct cloaking devices with varying bathymetry and a set of rigid vertical inclusions. Finally, we showed numerically that such devices lead to a cloaking effect for a structure placed in the coastal zone for
$kh\in [0.2;0.7]$
for linear water waves, with a reduction of the mean drift force. Although our study has been carried out in the MSE approximation, we hope that it will foster theoretical and numerical effort in water wave cloaking in the framework of the full linear theory, and that it will serve as a basis for experimental studies in linear and nonlinear water wave cloaking.
