1. Introduction
The method of transformation media offers an unprecedented way of steering a variety of waves through designing the propagation medium. The scheme was first proposed by Leonhardt (Reference Leonhardt2006) and Pendry, Schurig & Smith (Reference Pendry, Schurig and Smith2006) for electromagnetic waves, but since then has been implemented also for other types of waves, such as acoustic waves (Chen & Chan Reference Chen and Chan2007; Cummer & Schurig Reference Cummer and Schurig2007; Zhang, Xia & Fang Reference Zhang, Xia and Fang2011; Huang, Zhong & Liu Reference Huang, Zhong and Liu2014), seismic waves (Brûlé et al. Reference Brûlé, Javelaud, Enoch and Guenneau2014), elastodynamic waves (Farhat et al. Reference Farhat, Guenneau, Enoch and Movchan2009), matter waves (Zhang et al. Reference Zhang, Genov, Sun and Zhang2008) and water waves (Chen et al. Reference Chen, Yang, Zi and Chan2009; Berraquero et al. Reference Berraquero, Maurel, Petitjeans and Pagneux2013). The method of transformation media establishes that, in any wave system described by a form-invariant governing equation, the trajectory of wave rays can be produced along any desirable path by changing the properties of the propagation medium. Depending on the complexity of the desired trajectory, the required change in the properties of the propagation medium may be quite involved. In fact, material properties found using the transformation media method are usually anisotropic and heterogeneous in space (Pendry et al. Reference Pendry, Schurig and Smith2006), and as a result may need to be artificially engineered (such a medium is sometimes called ‘metamaterial’ (Schurig et al. Reference Schurig, Mock, Justice, Cummer, Pendry, Starr and Smith2006)).
The design procedure based on the method of transformation media first considers a plain original space
$S_{o}$
, with constant medium properties, in which wave rays propagate along straight lines. Then a spatial coordinate transformation is found that transforms these trajectories (i.e. straight wave rays in
$S_{o}$
) to the desired (potentially curved and complex) trajectories which are expected in the physical space
$S_{p}$
. Note that, while the transformed governing equation guarantees that the wave rays travel along those desired paths, the equation itself does not merit any physical wave system. Nevertheless, since the governing equation of the wave system is form-invariant, the transformed equation must have the same form as the one in the original space, but with (potentially) different coefficients. Noting that the coefficients of the governing equation in the original space
$S_{o}$
are medium properties, it can be shown that the new coefficients in the transformed equation may be absorbed into these medium properties. Therefore, the transformed governing equation can be viewed as the governing equation in physical space describing wave propagation in a medium whose properties are different from those in the original space. Basically, the transformation is absorbed in the medium properties. Waves propagating in such a medium travel along those specific desired trajectories.
The most salient application of the method of transformation media is the design of invisibility cloaks. A cloak is a two- or three-dimensional patch that encloses the object which is to be made invisible. A perfect cloak must not reflect nor scatter any part of the incident wave, and must perfectly reconstruct the incident wave downstream of the object in such a way that waves outside the cloak bear no information of the presence of the cloaking patch nor the encased object. A cloak is also required to be omnidirectional, that is, to cloak itself and the enclosed object from incident waves independent of their direction of incidence (therefore, two periscopes installed back to back are not considered to form a cloak).
Since Maxwell’s equations are form-invariant, a cloak of invisibility for objects from electromagnetic waves can, in theory, be designed using the transformation media method. Such a cloak requires spatial variability of those material properties that appear in the Maxwell’s equations, i.e. permittivity and permeability (Leonhardt Reference Leonhardt2006; Pendry et al. Reference Pendry, Schurig and Smith2006). The governing equations for linear shallow-water waves have the same form as the single polarization Maxwell’s equation. In fact, there is a one-to-one mapping between the two equations and boundary conditions at the interface of two media. Therefore, in analogy to the electromagnetic waves, the idea of transformation media can be implemented for shallow-water waves. The issue, however, is that for shallow-water waves the two properties that need to be spatially variable are the water depth and the gravitational acceleration. The latter clearly cannot be altered.
Here, we introduce a nonlinear transformation via which the required spatial variability of the gravitational acceleration is transformed to the water depth variations. Through this nonlinear transformation, an (asymptotically) perfect and physically realizable cloak of invisibility is designed. The proposed scheme rigorously satisfies the required conditions on the governing equation, and we show that, if the cloak is large enough, the desired boundary conditions are also obtained. We elaborate details of the method of transformation media and design procedure in general, and the design of an invisibility cloak for shallow-water waves in particular. The designed cloak works for all frequency waves that satisfy the shallow-water wave assumption. We also present direct simulation of the cloak using the finite element method.
We would like to comment that cloaking of surface gravity waves has also been investigated via alternative methodologies. Farhat et al. (Reference Farhat, Enoch, Guenneau and Movchan2008) use concentric arrays of surface-piercing cylinders to cloak an area from surface gravity–capillary waves. Through homogenization theory, they show that the effect of the surrounding cylinders appears as two different shear viscosities in two horizontal directions. These viscosities make the waves faster in one direction than the other, via which bending around the central area is made possible. In order to have a perfect cloak, however, besides a variable apparent viscosity, a radially variable density is also required that cannot be achieved physically but is approximated in Farhat et al. (Reference Farhat, Enoch, Guenneau and Movchan2008) by modifying the effective shear viscosity. This modification (also known as the reduced model) changes both the governing equation and the boundary conditions in such a way that none of them are fully satisfied. Such an approximate cloak nevertheless shows a relatively good performance, as appears from the observed surface profile (Farhat et al. Reference Farhat, Enoch, Guenneau and Movchan2008).
In Newman (Reference Newman2014), cloaking is investigated in deep water by surrounding a central cylinder (to be cloaked) by an array of surface-piercing (smaller) cylinders and also by a continuous ring. The geometries of the surrounding bodies are then optimized towards a cloak by minimizing the scattered energy. Calculation of the scattering pattern and the optimizations were performed by commercial software. The cloak presented by Newman (Reference Newman2014) is non-axisymmetric and hence is not omnidirectional. In fact Newman (Reference Newman2014) states that a perfect cloaking can only be achieved with non-axisymmetric structures. In a similar effort, Porter & Newman (Reference Porter and Newman2014) define a topographic ring about a central cylinder, and by minimizing the scattered wave energy converge towards a cloak. Likewise in this cloak the highest efficiency cloaking is obtained for non-axisymmetric topographies, which makes the cloak directional.
Alam (Reference Alam2012) achieved cloaking of surface objects in two-layer density-stratified fluids using seabed ripple invoking Bragg resonance between surface and interfacial waves (cf. Alam, Liu & Yue Reference Alam, Liu and Yue2009a ,Reference Alam, Liu and Yue b ). In this approach, waves are transported, upstream of a floating object, from the surface to the interface where waves pass underneath the surface object. On the downstream side of the object, interfacial waves can then be moved back to the surface through a similar mechanism. Downstream waves bypass the surface object and hence bear no trace of it.
2. Theory of the transformation media method
The governing equations for the propagation of long waves (i.e.
$kh\ll 1$
, where
$k$
is the wavenumber and
$h=h(\boldsymbol{x})$
is the water depth) in an inviscid, irrotational and incompressible fluid, i.e. potential flow, are
$$\begin{eqnarray}\displaystyle & \partial \boldsymbol{u}/\partial t+g\boldsymbol{{\rm\nabla}}{\it\eta}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \partial {\it\eta}/\partial t+\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }(h\boldsymbol{u})=0, & \displaystyle\end{eqnarray}$$
$\boldsymbol{u}$
represents the velocity vector,
${\it\eta}(\boldsymbol{x},t)$
is the water surface elevation measured from the mean water level,
$g$
is the gravitational acceleration, and
$\boldsymbol{{\rm\nabla}}$
is the horizontal gradient operator. Equations (2.1a
) and (2.1b
) are, respectively, the momentum and mass balance equations. Combining these two equations, assuming waves are time-harmonic with frequency
${\it\omega}$
and gravitational acceleration
$g$
is constant, we obtain the well-known shallow-water wave equation (e.g. Mei, Stiassnie & Yue Reference Mei, Stiassnie and Yue2005) 
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\left(h\boldsymbol{{\rm\nabla}}{\it\eta}\right)+\frac{{\it\omega}^{2}}{g}{\it\eta}=0.\end{eqnarray}$$
Equation (2.2) is form-invariant, i.e. it keeps its form under any arbitrary coordinate transformation. To show this property, we integrate (2.2) over a specific (yet arbitrary) volume
${\it\Omega}$
whose surface is denoted by
$\partial {\it\Omega}$
, and use the divergence theorem to get
$$\begin{eqnarray}\int _{\partial {\it\Omega}}h\boldsymbol{{\rm\nabla}}{\it\eta}\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S+\int _{{\it\Omega}}\frac{{\it\omega}^{2}}{g}{\it\eta}\,\text{d}V=0,\end{eqnarray}$$
where
$\boldsymbol{n}=(n_{1},n_{2})$
is the normal vector to the surface
$\partial {\it\Omega}$
, and
$\text{d}S$
and
$\text{d}V$
are, respectively, infinitesimal surface and volume elements. We now consider a coordinate transformation in the form
$\boldsymbol{x}^{\prime }=\boldsymbol{x}^{\prime }(\boldsymbol{x})$
with the transformation Jacobian
$\unicode[STIX]{x1D641}=\{\partial x_{i}^{\prime }/\partial x_{j}\}$
, where primed variables are in the new coordinate system. Under this transformation, (2.3) turns into
$$\begin{eqnarray}\int _{\partial {\it\Omega}^{\prime }}h\frac{\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}}{|\unicode[STIX]{x1D641}|}\boldsymbol{{\rm\nabla}}^{\prime }{\it\eta}\boldsymbol{\cdot }\boldsymbol{n}^{\prime }\,\text{d}S^{\prime }+\int _{{\it\Omega}^{\prime }}\frac{1}{|\unicode[STIX]{x1D641}|}\,\frac{{\it\omega}^{2}}{g}{\it\eta}\,\text{d}V^{\prime }=0,\end{eqnarray}$$
where we have used the relations
$\text{d}V^{\prime }=|\unicode[STIX]{x1D641}|\,\text{d}V$
and
$\boldsymbol{n}^{\prime }\,\text{d}S^{\prime }=|\unicode[STIX]{x1D641}|\unicode[STIX]{x1D641}^{-\text{T}}\boldsymbol{n}\,\text{d}S$
in which
$|\unicode[STIX]{x1D653}|$
stands for the determinant of the matrix
$\unicode[STIX]{x1D653}$
. With the help of the divergence theorem on (2.4), and noting that this equation is valid for any arbitrary volume, we obtain
$$\begin{eqnarray}\boldsymbol{{\rm\nabla}}^{\prime }\boldsymbol{\cdot }\left(\unicode[STIX]{x1D65D}^{\prime }\boldsymbol{{\rm\nabla}}^{\prime }{\it\eta}\right)+\frac{{\it\omega}^{2}}{g^{\prime }}{\it\eta}=0,\end{eqnarray}$$
where
$\unicode[STIX]{x1D65D}^{\prime }=h\,\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}/|\unicode[STIX]{x1D641}|$
and
$g^{\prime }=|\unicode[STIX]{x1D641}|g$
. Clearly (2.5) in the primed coordinates has the same form as (2.2) in the original coordinates, with the apparent values of
$h,g$
in the original coordinates replaced by
$\unicode[STIX]{x1D65D}^{\prime },g^{\prime }$
in the transformed coordinates. Therefore the governing equation (2.2) for long waves is form-invariant under coordinate transformation. This property can also be shown by a one-to-one mapping between the long-wave equation and the two-dimensional Maxwell’s equation for a transverse magnetic (TM) polarization wave in an isotropic medium with constant permittivity.
With the knowledge that (2.2) is form-invariant, we can invoke the transformation media method to steer long gravity waves in any desired direction. At this stage it is helpful to use some notation from differential geometry to express our form-invariant governing equation (2.2). It can be shown that for any coordinate system (
$x^{1},x^{2}$
) with the metric tensor
$\unicode[STIX]{x1D642}$
, (2.2) can be rewritten as (see e.g. Young Reference Young1992)
$$\begin{eqnarray}\frac{1}{\sqrt{|\unicode[STIX]{x1D642}|}}\,\boldsymbol{\partial }\boldsymbol{\cdot }\left(\sqrt{|\unicode[STIX]{x1D642}|}\,h\,\unicode[STIX]{x1D642}^{-1}\,\boldsymbol{\partial }{\it\eta}\right)+\frac{{\it\omega}^{2}}{g}{\it\eta}=0,\end{eqnarray}$$
where
$\boldsymbol{\partial }=(\partial /\partial x^{1},\partial /\partial x^{2})$
is the partial derivative operator. For instance, if the coordinate system of interest is Cartesian, its metric tensor
$\unicode[STIX]{x1D642}_{c}=\unicode[STIX]{x1D644}$
is an identity matrix and the traditional form of (2.2) is readily obtained.
We first consider an original space
$S_{o}$
in which wave rays are along specific trajectories. For example, if the water depth is constant and we have long-crested waves, then wave rays are along straight parallel lines. We define, in this space, a coordinate system
$(x_{o}^{1},x_{o}^{2})$
with the metric tensor
$\unicode[STIX]{x1D642}_{o}$
. According to (2.6), equation (2.2) in this space is written as
$$\begin{eqnarray}\frac{1}{\sqrt{|\unicode[STIX]{x1D642}_{o}|}}\,\boldsymbol{\partial }_{o}\boldsymbol{\cdot }\left(\sqrt{|\unicode[STIX]{x1D642}_{o}|}\,h\,\unicode[STIX]{x1D642}_{o}^{-1}\,\boldsymbol{\partial }_{o}{\it\eta}\right)+\frac{{\it\omega}^{2}}{g}{\it\eta}=0,\end{eqnarray}$$
with
$\boldsymbol{\partial }_{o}=(\partial /\partial x_{o}^{1},\partial /\partial x_{o}^{2})$
.
We then find a transformation
$T$
that maps the original space
$S_{o}$
to the physical space
$S_{p}$
in such a way that the wave ray trajectories in
$S_{o}$
are transformed to the desired trajectories in
$S_{p}$
. This transformation transforms the coordinate
$(x_{o}^{1},x_{o}^{2})$
to a new coordinate
$(x_{p}^{1},x_{p}^{2})$
whose metric is
$\unicode[STIX]{x1D642}_{p}$
and as a result (2.7) becomes
$$\begin{eqnarray}\frac{1}{\sqrt{|\boldsymbol{Q}|}}\,\boldsymbol{\partial }_{p}\boldsymbol{\cdot }\left(\sqrt{|\boldsymbol{Q}|}\,h\,\boldsymbol{Q}^{-1}\,\boldsymbol{\partial }_{p}{\it\eta}\right)+\frac{{\it\omega}^{2}}{g}{\it\eta}=0,\end{eqnarray}$$
where
$\boldsymbol{\partial }_{p}=(\partial /\partial x_{p}^{1},\partial /\partial x_{p}^{2})$
and
$\boldsymbol{Q}^{-1}=\unicode[STIX]{x1D641}\unicode[STIX]{x1D642}_{o}^{-1}\unicode[STIX]{x1D641}^{\text{T}}$
in which
$\unicode[STIX]{x1D641}\equiv \{\partial x_{p}^{i}/\partial x_{o}^{j}\}$
is the Jacobian of the transformation
$T$
. Note that in (2.8)
$\boldsymbol{Q}\neq \unicode[STIX]{x1D642}_{p}$
, i.e.
$\boldsymbol{Q}$
is not the metric tensor of the coordinate system
$(x_{p}^{1},x_{p}^{2})$
and therefore wave rays under (2.8) move along different trajectories than those described by (2.7). Basically, (2.8) is not the governing equation for water waves in a physical space.
Now consider the propagation of long waves in water of depth
$\tilde{\unicode[STIX]{x1D65D}}$
and under the action of gravitational acceleration
$\tilde{g}$
in the physical space, i.e.
$$\begin{eqnarray}\frac{1}{\sqrt{|\unicode[STIX]{x1D642}_{p}|}}\,\boldsymbol{\partial }_{p}\boldsymbol{\cdot }\left(\sqrt{|\unicode[STIX]{x1D642}_{p}|}\,\tilde{\unicode[STIX]{x1D65D}}\,\unicode[STIX]{x1D642}_{p}^{-1}\,\boldsymbol{\partial }_{p}{\it\eta}\right)+\frac{{\it\omega}^{2}}{\tilde{g}}{\it\eta}=0.\end{eqnarray}$$
Nevertheless, if we choose specific values for
$\tilde{\unicode[STIX]{x1D65D}}$
and
$\tilde{g}$
as
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D65D}}=\frac{\sqrt{|\unicode[STIX]{x1D642}_{o}|}}{\sqrt{|\unicode[STIX]{x1D642}_{p}|}}\frac{\unicode[STIX]{x1D641}\unicode[STIX]{x1D642}_{o}^{-1}\unicode[STIX]{x1D641}^{\text{T}}\unicode[STIX]{x1D642}_{p}}{|\unicode[STIX]{x1D641}|}h,\quad \tilde{g}=\frac{\sqrt{|\unicode[STIX]{x1D642}_{p}|}\,|\unicode[STIX]{x1D641}|}{\sqrt{|\unicode[STIX]{x1D642}_{o}|}}g,\end{eqnarray}$$
then (2.9) and (2.8) become identical. Physically speaking, this means that if in the physical space the water depth and the gravitational acceleration are given by (2.10), then the incident wave will be travelling along desired trajectories given by the transformation
$T$
. Note that (2.10) in a Cartesian coordinate system simplifies to
$\tilde{\unicode[STIX]{x1D65D}}=h\,\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}/|\unicode[STIX]{x1D641}|$
and
$\tilde{g}=|\unicode[STIX]{x1D641}|g$
(see also appendix A).
3. Implementation for shallow-water waves
There are two major issues in realizing conditions given by (2.10) in the physical world: (i) in (2.10) the water depth
$\tilde{\unicode[STIX]{x1D65D}}$
is a tensor while water depth in the real world is a scalar, and (ii) the gravitational acceleration
$\tilde{g}$
as is required by (2.10) is different from the physical world’s gravitational acceleration
$g$
.
Generally speaking in (2.2) water depth
$h$
can be regarded as a tensor. In the simplest case of a constant water depth
$h(x,y)=h_{0}$
, the tensor elements are
$\unicode[STIX]{x1D629}_{11}=\unicode[STIX]{x1D629}_{22}=h_{0}$
and
$\unicode[STIX]{x1D629}_{12}=\unicode[STIX]{x1D629}_{21}=0$
. However, the challenge of issue (i) is that in (2.10) the water depth tensor may have non-identical main diagonal elements. For example, if the physical space
$S_{p}$
is described by a Cartesian coordinate system, then the required water depth
$\tilde{\unicode[STIX]{x1D65D}}$
in (2.10) is given by
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D65D}_{c}}=\sqrt{|\unicode[STIX]{x1D642}_{o}|}\frac{\unicode[STIX]{x1D641}\unicode[STIX]{x1D642}_{o}^{-1}\unicode[STIX]{x1D641}^{\text{T}}}{|\unicode[STIX]{x1D641}|}h,\end{eqnarray}$$
which is a symmetric tensor whose elements depend on
$\unicode[STIX]{x1D641}$
and
$\unicode[STIX]{x1D642}_{o}$
.
To address this issue, consider a spatially varying one-dimensional topography
$h=h(x)$
with the characteristic wavelength
$L_{x}\ll {\it\lambda}$
, where
${\it\lambda}$
is the wavelength of the incident wave. Homogenization theory implemented on this problem yields an equivalent depth
$h_{x}$
in the
$x$
direction and an equivalent
$h_{y}$
in the
$y$
direction (e.g. Jikov, Oleinik & Kozlov Reference Jikov, Oleinik and Kozlov1994), or equivalently the tensor elements are
$\unicode[STIX]{x1D629}_{11}=h_{x},\unicode[STIX]{x1D629}_{22}=h_{y}$
and
$\unicode[STIX]{x1D629}_{12}=\unicode[STIX]{x1D629}_{21}=0$
.
Now note that in (3.1) the water depth tensor
$\tilde{\unicode[STIX]{x1D65D}}_{c}$
is symmetric and therefore it has two orthogonal eigenvectors of
$\boldsymbol{v}_{1}$
and
$\boldsymbol{v}_{2}$
with eigenvalues of
$v_{1}$
and
$v_{2}$
, respectively. In the basis of
$\boldsymbol{v}_{1}$
and
$\boldsymbol{v}_{2}$
, the water depth tensor
$\tilde{\unicode[STIX]{x1D65D}}_{c}$
is diagonal with elements
$v_{1}$
and
$v_{2}$
. If we define our
$x$
and
$y$
along
$\boldsymbol{v}_{1}$
and
$\boldsymbol{v}_{2}$
, then by a proper choice of
$h=h(x)$
, we can achieve
$\unicode[STIX]{x1D629}_{11}=v_{1}$
and
$\unicode[STIX]{x1D629}_{22}=v_{2}$
and therefore the required water depth tensor
$\tilde{\unicode[STIX]{x1D65D}}_{c}$
is obtained.
If the physical space
$S_{p}$
is described by any arbitrary (potentially curvilinear) coordinate system, related to the Cartesian coordinate by the mapping
${\it\bf\Lambda}$
, then the water depth in the new system
$\tilde{\unicode[STIX]{x1D65D}}_{new}$
is described by
$\tilde{\unicode[STIX]{x1D65D}}_{new}={\it\bf\Lambda}^{-1}\tilde{\unicode[STIX]{x1D65D}}_{c}{\it\bf\Lambda}$
. Clearly
$\tilde{\unicode[STIX]{x1D65D}}_{new}$
has eigenvectors of
${\it\bf\Lambda}^{-1}\boldsymbol{v}_{1}$
and
${\it\bf\Lambda}^{-1}\boldsymbol{v}_{2}$
with eigenvalues of
$v_{1}$
and
$v_{2}$
. In other words, the eigenvectors of the water depth in the new coordinate system are transformation of the eigenvectors of the Cartesian coordinate system, but the eigenvalues of the two systems are the same. Therefore, if
$\tilde{\unicode[STIX]{x1D65D}}_{new}$
is given, then
$h(x)$
can be determined by the eigenvalues of
$\tilde{\unicode[STIX]{x1D65D}}_{new}$
, with the direction of
$x$
and
$y$
being along the eigenvectors of
$\tilde{\unicode[STIX]{x1D65D}}_{new}$
.
The second issue has to do with the fact that, in (2.10), the gravitational acceleration
$\tilde{g}$
is not the same as the physical world’s gravitational acceleration
$g$
. Clearly, if
$\sqrt{|\unicode[STIX]{x1D642}_{p}|}\,|\unicode[STIX]{x1D641}|/\sqrt{|\unicode[STIX]{x1D642}_{o}|}=1$
, then this issue is resolved. However, this cannot always be achieved depending on the transformation required. Nevertheless, if the coefficient of
$g$
is (2.10) is any arbitrary constant, i.e.
$\sqrt{|\unicode[STIX]{x1D642}_{p}|}\,|\unicode[STIX]{x1D641}|/\sqrt{|\unicode[STIX]{x1D642}_{o}|}=C$
, then both terms of (2.9) can be multiplied by
$C$
, and then
$C$
can be absorbed into
$\tilde{\unicode[STIX]{x1D65D}}$
. Therefore, by defining a new water depth
$\bar{\unicode[STIX]{x1D65D}}=C\tilde{\unicode[STIX]{x1D65D}}$
and
$\bar{g}=g$
, we retrieve the required equation in the physical space. The significance of the arbitrary constant
$C$
is that it provides an additional degree of freedom which, as will be shown later, is critical in designing a specific set of transformations. It is to be noted that multiplying the equation by a constant, although it keeps the governing equation intact, affects the required boundary conditions. This is elaborated in detail in § 5.
In cases where the coefficient of
$g$
in (2.9) is neither unity nor a constant, an approximate approach is to use the following water depth
$\check{\unicode[STIX]{x1D65D}}$
and the gravitational acceleration
${\check{g}}$
:
$$\begin{eqnarray}\check{\unicode[STIX]{x1D65D}}=\unicode[STIX]{x1D641}\unicode[STIX]{x1D642}_{o}^{-1}\unicode[STIX]{x1D641}^{\text{T}}\unicode[STIX]{x1D642}_{p}\,h,\quad {\check{g}}=g,\end{eqnarray}$$
via which the speed of long waves remains the same, i.e.
$\tilde{\unicode[STIX]{x1D65D}}\tilde{g}=\check{\unicode[STIX]{x1D65D}}{\check{g}}$
. The origin of this simplification is that the refraction angle of the wave rays in oblique incidence at a depth discontinuity with constant gravity is a function of the ratio of the long-wave speed (i.e.
$\sqrt{gh}$
) in the two media (Mei et al.
Reference Mei, Stiassnie and Yue2005). But note that if the gravitational acceleration is spatially varying, which is the case here, the governing equation is not fully satisfied. This is a typical approximation widely used in the electromagnetic community to keep the permeability (or permittivity) constant (Cummer et al.
Reference Cummer, Popa, Schurig, Smith and Pendry2006; Schurig et al.
Reference Schurig, Mock, Justice, Cummer, Pendry, Starr and Smith2006; Cai et al.
Reference Cai, Chettiar, Kildishev and Shalaev2007) and the approximate nature of that has been noted (Cummer et al.
Reference Cummer, Popa, Schurig, Smith and Pendry2006; Yan, Ruan & Qiu Reference Yan, Ruan and Qiu2007).
4. Cloak of invisibility
An invisibility cloak for water waves is an area of (potentially) variable topography that bends the wave rays about a specific target area that is to be cloaked. The cloaked area and objects therein do not scatter incident waves, and therefore are invisible to far-field observers.
Here we aim at the cloaking of a circular area
$A_{c}$
with radius
$a$
. One way to achieve this is to transform a circular region of radius
$b>a$
(centred at the centre of
$A_{c}$
) to an annular ring of inner radius
$a$
and outer radius
$b$
centred again at the centre of
$A_{c}$
. Basically, the transformation must satisfy the conditions
$r_{p}(r_{o}=0)=a$
and
$r_{p}(r_{o}=b)=b$
, where subscripts
$o$
and
$p$
refer to the original and physical spaces, respectively.
One family of such transformations is
$$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}r_{p}=f(r_{o}),\quad 0\leqslant r_{o}\leqslant b,\quad \\ {\it\theta}_{p}={\it\theta}_{o},\quad \end{array}\right.\end{eqnarray}$$
where
$f(0)=a$
and
$f(b)=b$
. The region beyond the circle
$b$
(i.e.
$r_{o}>b$
) is identically mapped from
$S_{o}$
to
$S_{p}$
. The Jacobian of the transformation (4.1) is
$$\begin{eqnarray}\unicode[STIX]{x1D641}=\frac{\partial \left(r_{p},{\it\theta}_{p}\right)}{\partial \left(r_{o},{\it\theta}_{o}\right)}=\left(\begin{array}{@{}cc@{}}f^{\prime } & 0\\ 0 & 1\end{array}\right),\end{eqnarray}$$
where
$f^{\prime }=\text{d}f(r_{o})/\text{d}r_{o}$
. For the problem at hand, it is more convenient to use a polar coordinate system in both the original and physical spaces, with the following metric tensors
$$\begin{eqnarray}\unicode[STIX]{x1D642}_{o}=\left(\begin{array}{@{}cc@{}}1 & 0\\ 0 & r_{o}^{2}\end{array}\right),\quad \unicode[STIX]{x1D642}_{p}=\left(\begin{array}{@{}cc@{}}1 & 0\\ 0 & r_{p}^{2}\end{array}\right).\end{eqnarray}$$
Substituting (4.2) and (4.3) into (2.10), we obtain
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D65D}}=\left(\begin{array}{@{}cc@{}}\displaystyle \frac{f^{\prime }r_{o}}{r_{p}} & 0\\ 0 & \displaystyle \frac{r_{p}}{f^{\prime }r_{o}}\end{array}\right)h,\quad \tilde{g}=\frac{f^{\prime }r_{p}}{r_{o}}\,g.\end{eqnarray}$$
Equation (4.4) gives the values of
$\tilde{\unicode[STIX]{x1D65D}}$
and
$\tilde{g}$
in the physical space in such a way that the wave rays follow the path determined by the transformation in (4.1). In order to find a mathematically perfect as well as a physically realizable cloak, we have to fulfill an additional condition that the ratio
$f^{\prime }r_{p}/r_{o}$
(i.e. coefficient of
$g$
in (4.4)) is constant. If this condition is satisfied, then as discussed in § 3, a new water depth
$\bar{\unicode[STIX]{x1D65D}}$
and physically admissible
$\bar{g}=g$
are found that steer the wave rays along the trajectories defined by the transformation (4.1). For the considered case of a circular cloak it is easy to show that the following nonlinear transformation
$$\begin{eqnarray}\left\{\begin{array}{@{}l@{}}r_{p}=f(r_{o})=\sqrt{\left(1-\displaystyle \frac{a^{2}}{b^{2}}\right)r_{o}^{2}+a^{2}},\quad 0\leqslant r_{o}\leqslant b,\quad \\ {\it\theta}_{p}={\it\theta}_{o},\quad \end{array}\right.\end{eqnarray}$$
satisfies the aforementioned condition and gives rise to
$$\begin{eqnarray}\tilde{\unicode[STIX]{x1D65D}}=\left(\begin{array}{@{}cc@{}}\displaystyle \frac{r_{p}^{2}-a^{2}}{r_{p}^{2}} & 0\\ 0 & \displaystyle \frac{r_{p}^{2}}{r_{p}^{2}-a^{2}}\end{array}\right)h,\quad \tilde{g}=\left(1-\frac{a^{2}}{b^{2}}\right)\,g,\end{eqnarray}$$
and hence the final form of the water depth and gravitational acceleration in the physical space is obtained as
$$\begin{eqnarray}\bar{\unicode[STIX]{x1D65D}}=\left(1-\frac{a^{2}}{b^{2}}\right)\left(\begin{array}{@{}cc@{}}\displaystyle \frac{r_{p}^{2}-a^{2}}{r_{p}^{2}} & 0\\ 0 & \displaystyle \frac{r_{p}^{2}}{r_{p}^{2}-a^{2}}\end{array}\right)h,\quad \bar{g}=\,g.\end{eqnarray}$$
Water waves propagating over a topography given by (4.7) in the physical world of gravitational acceleration
$g$
move along trajectories given by (4.1) which are also shown in figure 1. Clearly wave rays go around the cylinder in the middle and hence do not interact with that cylinder. The cylinder is now cloaked. Note that although (4.7) fully satisfies the transformed equation, it changes the flux boundary condition at the cloak’s boundary by a factor of
$(1-a^{2}/b^{2})$
. Clearly, as
$a/b\rightarrow 0$
, i.e. as the cloak size increases, the desired boundary condition is asymptotically satisfied. This is discussed in more detail in § 5 where we derive the analytical solution of this cloak.
Trajectories of the wave rays due to transformation in (4.1). These trajectories show the flow path of the wave’s energy.
In the special case when both the original and physical spaces are described by a Cartesian coordinate system, which is typically the case in the cloaking literature, the water depth tensor and gravitational acceleration (4.6) can be obtained directly, as elaborated in appendix A. We would also like to note that ‘rays’ and ‘trajectories’ are used here merely to discuss the physical picture, and do not imply the use of the ‘ray theory’.
5. Analytical solution
In this section, we present an analytical solution to the governing equation (2.2) for the incoming wave and the scattered field as waves travel over the topography given by (4.7). Equation (2.2) in polar coordinates reads
$$\begin{eqnarray}\frac{1}{r}\frac{\partial }{\partial r}\left(r\unicode[STIX]{x1D629}_{r}\frac{\partial {\it\eta}}{\partial r}\right)+\frac{1}{r^{2}}\frac{\partial }{\partial {\it\theta}}\left(\unicode[STIX]{x1D629}_{{\it\theta}}\frac{\partial {\it\eta}}{\partial {\it\theta}}\right)+\frac{{\it\omega}^{2}}{g}{\it\eta}=0,\end{eqnarray}$$
where we have assumed the water depth tensor to be diagonal, which is the case for the cylindrical cloaks considered here. Assuming that the solution to (5.1) is separable, i.e.
${\it\eta}=\text{Re}\{{\it\Psi}(r){\it\Theta}({\it\theta})\}$
, we obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}^{2}}{\text{d}{\it\theta}^{2}}{\it\Theta}+m^{2}{\it\Theta}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}r}\left(r\unicode[STIX]{x1D629}_{r}\frac{\text{d}{\it\Psi}}{\text{d}r}\right)+\frac{{\it\omega}^{2}r}{g}{\it\Psi}-\frac{m^{2}}{r}\unicode[STIX]{x1D629}_{{\it\theta}}{\it\Psi}=0. & \displaystyle\end{eqnarray}$$
${\it\Theta}({\it\theta})=\exp (\text{i}m{\it\theta})$
. Since
$\unicode[STIX]{x1D629}_{r}$
and
$\unicode[STIX]{x1D629}_{{\it\theta}}$
are different in the two regions
$a\leqslant r\leqslant b$
and
$r\geqslant b$
, we break (5.3) into two equations, each for one of the domains, 
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}r}\left(\frac{r^{2}-a^{2}}{r}h\frac{\text{d}{\it\Psi}}{\text{d}r}\right)+\frac{{\it\omega}^{2}br}{g(b^{2}-a^{2})}{\it\Psi}-\frac{m^{2}r}{r^{2}-a^{2}}h{\it\Psi}=0,\quad a\leqslant r\leqslant b, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}}{\text{d}r}\left(rh\frac{\text{d}{\it\Psi}}{\text{d}r}\right)+\frac{{\it\omega}^{2}r}{g}{\it\Psi}-\frac{m^{2}}{r}h{\it\Psi}=0,\quad r\geqslant b, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}{\it\Psi}(r)=\left\{\begin{array}{@{}ll@{}}{\it\beta}_{m}^{i}\text{J}_{m}(\tilde{k}\sqrt{r^{2}-a^{2}})+{\it\beta}_{m}^{s}\text{H}_{m}^{(1)}(\tilde{k}\sqrt{r^{2}-a^{2}}),\quad & a\leqslant r\leqslant b,\\ {\it\alpha}_{m}^{i}\text{J}_{m}(kr)+{\it\alpha}_{m}^{s}\text{H}_{m}^{(1)}(kr),\quad & r\geqslant b,\end{array}\right.\end{eqnarray}$$
where
$\text{J}_{m}$
and
$\text{H}_{m}^{(1)}$
are respectively order-
$m$
Bessel function and Hankel function of the first kind and
$k^{2}={\it\omega}^{2}/(gh)$
and
$\tilde{k}=kb/(b^{2}-a^{2})^{1/2}$
. The Hankel function of the first kind asymptotically approaches an outgoing wave as
$r$
tends to infinity, and represents the scattering field.
Noting that an incident plane wave of the form
${\it\eta}_{inc}=\text{Re}\{{\it\eta}_{0}\exp (\text{i}kx)\}$
can be expressed as
${\it\eta}_{inc}=\text{Re}\{{\it\eta}_{0}\sum _{m=-\infty }^{\infty }\text{i}^{m}\text{J}_{m}(kr)\text{e}^{\text{i}m{\it\theta}}\}$
, the coefficients
${\it\alpha}_{m}^{i}$
in (5.6) are readily obtained. The boundary condition on the surface of the cylinder
$r=a$
is the no-flux (Neumann) condition and implies
${\it\beta}_{m}^{s}=0$
. The other three coefficients are in general non-zero. The free surface elevation
${\it\eta}(r,{\it\theta})$
is obtained as
$$\begin{eqnarray}{\it\eta}(r,{\it\theta})=\text{Re}\left\{\begin{array}{@{}ll@{}}\displaystyle \mathop{\sum }_{m=-\infty }^{\infty }{\it\beta}_{m}^{i}\text{J}_{m}(\tilde{k}\sqrt{r^{2}-a^{2}})\text{e}^{\text{i}m{\it\theta}},\quad & a\leqslant r\leqslant b,\\ \displaystyle \mathop{\sum }_{m=-\infty }^{\infty }({\it\alpha}_{m}^{i}\text{J}_{m}(kr)+{\it\alpha}_{m}^{s}\text{H}_{m}^{(1)}(kr))\text{e}^{\text{i}m{\it\theta}},\quad & r\geqslant b.\end{array}\right.\end{eqnarray}$$
At
$r=b$
, both surface elevation
${\it\eta}$
and mass flux
$\unicode[STIX]{x1D629}_{r}u_{r}=(\text{i}g\unicode[STIX]{x1D629}_{r}/{\it\omega})\partial {\it\eta}/\partial r$
must be continuous. Using the parameters of the cloak given by (4.7) we obtain from (5.6) the following two relations:
$$\begin{eqnarray}\displaystyle & {\it\alpha}_{m}^{i}\text{J}_{m}(kb)+{\it\alpha}_{m}^{s}\text{H}_{m}^{(1)}(kb)={\it\beta}_{m}^{i}\text{J}_{m}(kb), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\alpha}_{m}^{i}\text{J}_{m}^{\prime }(kb)+{\it\alpha}_{m}^{s}\text{H}_{m}^{\prime (1)}(kb)=\left(1-\frac{a^{2}}{b^{2}}\right){\it\beta}_{m}^{i}\text{J}_{m}^{\prime }(kb), & \displaystyle\end{eqnarray}$$
where the prime on
$\text{J}_{m}$
and
$\text{H}_{m}^{(1)}$
stands for the derivative with respect to their arguments. Therefore
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\beta}_{m}^{i}=\frac{\text{J}_{m}^{\prime }(kb)\text{H}_{m}^{(1)}(kb)-\text{J}_{m}(kb)\text{H}_{m}^{\prime (1)}(kb)}{\displaystyle \left(1-\frac{a^{2}}{b^{2}}\right)\text{J}_{m}^{\prime }(kb)\text{H}_{m}^{(1)}(kb)-\text{J}_{m}(kb)\text{H}_{m}^{\prime (1)}(kb)}{\it\alpha}_{m}^{i}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\alpha}_{m}^{s}=\frac{\displaystyle \frac{a^{2}}{b^{2}}\text{J}_{m}(kb)\text{J}_{m}^{\prime }(kb)}{\displaystyle \left(1-\frac{a^{2}}{b^{2}}\right)\text{J}_{m}^{\prime }(kb)\text{H}_{m}^{(1)}(kb)-\text{J}_{m}(kb)\text{H}_{m}^{\prime (1)}(kb)}{\it\alpha}_{m}^{i}. & \displaystyle\end{eqnarray}$$
The coefficient
${\it\alpha}_{m}^{s}$
that appears in front of the Hankel function in (5.6) shows the scattered field in the area
$r\geqslant b$
and therefore for a perfect cloak must be zero. It is easy to see that this is asymptotically obtained when
$b/a\rightarrow \infty$
. In other words, if
$b/a\gg 1$
then we asymptotically obtain
${\it\alpha}_{m}^{s}=0,{\it\beta}_{m}^{s}=0$
and
${\it\beta}_{m}^{i}={\it\alpha}_{m}^{i}$
(cf. figures 2 and 4).
Analytically obtained cloaking factor
$\mathscr{L}$
calculated from (5.7) and (5.14) as a function of
$b/a$
. The cloaking factor exponentially decays and asymptotically approaches the ideal cloaking factor of
$\mathscr{L}=0$
. At
$b/a=5$
, the cloaking factor
$\mathscr{L}<e^{-5}$
for both
${\it\lambda}/a=4$
and 8. In calculating
$\mathscr{L}$
, the series in (5.7) is truncated at
$m=\pm 40$
for which the value of
$\mathscr{L}$
is converged.
In order to quantitatively measure the efficiency of the implemented cloak, we calculate the scattering cross-section (measure of the total energy scattered to infinity) for the two cases studied here. The scattering cross-section
${\it\sigma}({\it\theta})$
can be defined as (e.g. Porter & Newman Reference Porter and Newman2014)
$$\begin{eqnarray}{\it\sigma}({\it\theta})=\lim _{r\rightarrow \infty }2{\rm\pi}r\frac{\overline{{\it\Phi}}_{sca}^{2}(r,{\it\theta})}{\overline{{\it\Phi}}_{inc}^{2}(r,{\it\theta})},\end{eqnarray}$$
where
${\it\Phi}_{sca}(r,{\it\theta})$
and
${\it\Phi}_{inc}(r,{\it\theta})$
are respectively the scattered wave velocity potential and the incident wave velocity potential. The velocity potential
${\it\Phi}$
is obtained from the dynamic boundary condition on the free surface, i.e.
${\it\Phi}=(\text{i}g/{\it\omega}){\it\eta}$
. The scattering field is calculated by subtracting the incident wave potential from the total velocity potential, i.e.
${\it\Phi}_{sca}={\it\Phi}_{tot}-{\it\Phi}_{inc}$
. The overbars stand for the r.m.s. values in time. The total scattering cross-section
${\it\Sigma}$
is therefore
$$\begin{eqnarray}{\it\Sigma}=\frac{1}{2{\rm\pi}}\oint {\it\sigma}({\it\theta})\,\text{d}{\it\theta}.\end{eqnarray}$$
The ratio of the total energy scattered by the object in the presence of the cloak
${\it\Sigma}_{clk}$
to the total energy scattered when the seabed is flat
${\it\Sigma}_{flt}$
is defined as the cloaking factor (Porter & Newman Reference Porter and Newman2014),
$$\begin{eqnarray}\mathscr{L}=\frac{{\it\Sigma}_{clk}}{{\it\Sigma}_{flt}}.\end{eqnarray}$$
Perfect cloaking is achieved for
$\mathscr{L}=0$
when no energy is scattered from the cylinder, and the cloaking factor for a flat seabed is unity.
With this definition of the cloaking factor
$\mathscr{L}$
, the effect of the ratio
$b/a$
on the cloaking efficiency can be quantitatively addressed. For instance, in figure 2 we present
$\mathscr{L}$
(for an enhanced accuracy averaged over the area
$8\leqslant r/a\leqslant 10$
) as a function of
$b/a$
and for two values of
${\it\lambda}/a=4$
and 8. As can be seen, with the increase in
$b/a$
, the cloaking factor
$\mathscr{L}$
exponentially decays to the ideal value of
$\mathscr{L}=0$
. For instance, for
$b/a=5$
,
$\mathscr{L}<\text{e}^{-5}$
for both
${\it\lambda}/a=4$
and 8.
6. Numerical results
For the implementation of the water depth tensor given by (4.7), as discussed in § 3, we employ homogenization theory on a spatially rapidly varying topography. To achieve this for the case of (4.7), consider a water depth
$h=h(r)$
that alternates in the radial direction between (constant) depths
$h_{1}$
and
$h_{2}$
with respective widths
$b_{1}$
and
$b_{2}$
. If
$b_{1},b_{2}\ll {\it\lambda}$
then homogenization theory gives
$$\begin{eqnarray}\unicode[STIX]{x1D629}_{r}^{eq}=\frac{h_{1}h_{2}(b_{1}+b_{2})}{b_{1}h_{2}+b_{2}h_{1}},\quad \unicode[STIX]{x1D629}_{{\it\theta}}^{eq}=\frac{b_{1}h_{1}+b_{2}h_{2}}{b_{1}+b_{2}}.\end{eqnarray}$$
By equating (6.1) with those depths required by (4.7),
$h_{1}(r_{p})$
and
$h_{2}(r_{p})$
are obtained. In our simulation, for simplicity, we consider
$b_{1}=b_{2}$
, for which we have
$h_{1,2}(r_{p})=\bar{\unicode[STIX]{x1D629}}_{22}[1\pm (1-\bar{\unicode[STIX]{x1D629}}_{11}/\bar{\unicode[STIX]{x1D629}}_{22})^{1/2}]$
.
We use finite element method with C++ library deal.II (Bangerth, Hartmann & Kanschat Reference Bangerth, Hartmann and Kanschat2007) for the direct simulation of the cloaking (for details on validation and convergence of the scheme, see e.g. Bangerth & Kanschat (Reference Bangerth and Kanschat1999), Bangerth et al. (Reference Bangerth, Hartmann and Kanschat2007), Carraro et al. (Reference Carraro, Goll, Marciniak-Czochra and Mikelić2013) and Riedlbauer et al. (Reference Riedlbauer, Drexler, Drummer, Steinmann and Mergheim2014)). We choose a domain of computation
$x,y\in [-4,4]$
. The mesh is generated using deal.II library with 2048 nodes in each direction. We chose
${\it\delta}t$
such that the Courant number in the simulations is
$C=0.5$
. We implement a numerical wavemaker on the left boundary. Other vertical boundaries are chosen as reflecting walls. The boundary condition on the inner cylinder is Neumann, which corresponds to no-flux condition.
The water surface profile is shown in figure 3(a–d) for two ratios of the wavelength to the cylinder radius
${\it\lambda}/a=4$
and 8 and for each case with and without the cloak. We have chosen maximum
$kh=0.2$
, which happens at the inner cylinder’s boundary, where we replaced the infinite depth with a depth
$O(10)$
times larger than the mean water depth outside the cloak. Specifically, in this paper, we have used 12 times larger than the mean depth. We would like to emphasize again that (2.2) is only valid in the limit of shallow water, i.e.
$kh\ll 1$
. For the investigated cases presented here, we have made sure that this condition is always satisfied throughout the domain of simulation, as described above.
Cloaking of a circular cylinder from surface gravity waves. The top and bottom panels show two cases of
${\it\lambda}/a=4$
and 8 respectively. (a,c) The water surface pattern (i.e. incident and scattered waves) in the vicinity of the cylinder if the seabed is flat. (b,d) The water surface with the implementation of the cloak (
$b/a=4$
). The outer perimeter of the cloak is marked by the black circle.
Figure 3(a,c) show the incident and scattered waves for respectively
${\it\lambda}/a=4$
and 8 and when the water depth is constant. Clearly a larger cylinder results in a higher scattering. The water surface profile with the implementation of the cloak of (4.7) with
$b/a=4$
is shown in figure 3(b,d). The black circle shows the outer radius of the annular cloak. As is seen from figure 3(b,d), the cloaking is visually perfect in both cases.
With the implementation of the cloak designed based on (4.7) with
$b/a=4$
, we obtained
$\mathscr{L}=0.013$
for
${\it\lambda}/a=4$
and
$\mathscr{L}=0.021$
for
${\it\lambda}/a=8$
. This means that
${\sim}99\,\%$
of the scattering energy is prevented by the proposed cloak. The exact theory of § 5 predicts
$\mathscr{L}=0.0078$
and
$\mathscr{L}=0.013$
for respectively
${\it\lambda}/a=4$
and 8. The discrepancy comes from computational errors, homogenization error, errors due to the finite time simulation, and the treatment of the singularity at
$r=a$
, which in numerical implementation is replaced by a large but a finite water depth chosen such that the shallow-water assumption is not violated. In the physical space, if the incident wavelength is
${\it\lambda}=200~\text{m}$
, the mean water depth outside the cloak is 50 cm and the maximum depth at
$r=a$
is taken as 6.3 m, the above two cases correspond to cylinder radii of
$a=50$
and
$25~\text{m}$
and outer cloak radii of
$b=200$
and
$100~\text{m}$
, respectively. If a cloaking factor of
$\mathscr{L}=0.1$
(
${\sim}90\,\%$
cloaking) is acceptable, then the same cloaking can be achieved over a mean water depth of 2 m.
In § 5, we derived analytically that the cloaking factor asymptotically approaches zero as the ratio
$b/a$
increases. This can also be investigated numerically. Figure 4 compares the water surface about a cylinder with a cloak implemented on the seabed for
$b/a=2$
(a) and
$b/a=4$
(b). While cloaking in figure 4(b) is almost perfect, the case of figure 4(a) is clearly not as good. The discrepancy, as discussed in § 5, comes from the discontinuity in the flux at the boundary
$r=b$
. This discontinuity exponentially fades away, and the cloaking efficiency enhances, as the ratio
$b/a$
increases.
Effect of the cloak size
$b/a$
on the cloaking efficiency: (a) for the cloak radius
$b/a=2$
for which
$\mathscr{L}=0.11$
and (b) (same as figure 3
b) for the cloak radius
$b/a=4$
for which
$\mathscr{L}=0.013$
. Clearly the case of
$b/a=4$
(b) achieves a much better cloaking efficiency than
$b/a=2$
(a). As discussed in § 5, cloaking efficiency, parametrized by the cloaking factor
$\mathscr{L}$
, exponentially improves as the ratio
$b/a$
increases.
Besides the specific form of the topography proposed and discussed here, many other forms of seabed variations may result in the desired depth tensor of (4.7). For instance, an anisotropic bed can be formed by overlaying two interlocking fine-meshed combs aligned in the radial and azimuthal directions which independently protrude to different height profiles. On this topography, waves travelling radially only see the height profile of the circular comb as they pass uninterrupted through the teeth of the radial comb whilst circular waves only see the height profile of the radial comb.
7. Conclusion
Here we present the design of a symmetric cloak for shallow-water waves via the scheme of transformation media. The method of transformation media states that if the governing equation of a wave system is form-invariant, then wave ray trajectories can be forced along any desired path by appropriately designing the properties of the medium of propagation. The challenge, for the case of water waves, is that usual linear implementations of the scheme require water depth and gravitational acceleration to be spatially variable along the path. The latter is simply impossible. Here we present a rigorous nonlinear alternative that relaxes the condition on the gravitational acceleration, and guarantees an asymptotically perfect cloak by mere variation in the seabed topography. As an example, a cloak for a circular cylinder is designed and its performance under incident waves is presented analytically and via a time-domain direct simulation. The presented cloak is asymptotically perfect, physically realizable, omnidirectional and works for all frequency waves that satisfy the shallow-water wave assumption.
Acknowledgements
This publication was made possible, in part, with the support from NSF (Grant No. CBET-1414579-EAGER), the American Bureau of Shipping, and the U.C. Berkeley Committee on Research. We would also like to thank the anonymous referees for valuable comments.
Appendix A
In § 2 we presented a general derivation of the method of transformation media that distinguishes between the Jacobian of a transformation
$\unicode[STIX]{x1D641}$
, and a metric tensor
$\unicode[STIX]{x1D642}$
. This derivation has the advantage that one can adopt any coordinate system for the physical space and/or the original space. For instance, in the problem of cylindrical cloaking, it is more convenient to use polar coordinates in both physical and original coordinate systems.
In the special case that we use Cartesian coordinates in both original and physical spaces, which is typically the common practice in the literature, the metric
$\unicode[STIX]{x1D642}$
is simply an identity tensor and (2.10) readily reduces to
$\tilde{\unicode[STIX]{x1D65D}}=h\,\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}/|\unicode[STIX]{x1D641}|$
and
$\tilde{g}=|\unicode[STIX]{x1D641}|g$
. For example, consider the transformation (4.1) that transforms the radial direction in physical and original space but keeps the angles intact. In a Cartesian coordinate system, this transformation can be written as
$$\begin{eqnarray}r=f(r^{\prime }),\quad r^{\prime }=\sqrt{x^{\prime 2}+y^{\prime 2}},\quad r=\sqrt{x^{2}+y^{2}},\end{eqnarray}$$
$$\begin{eqnarray}\cos {\it\theta}=\frac{x}{r}=\frac{x^{\prime }}{r^{\prime }}=\cos {\it\theta}^{\prime },\quad \sin {\it\theta}=\frac{y}{r}=\frac{y^{\prime }}{r^{\prime }}=\sin {\it\theta}^{\prime },\end{eqnarray}$$
whose Jacobian
$\unicode[STIX]{x1D641}$
is
$$\begin{eqnarray}\unicode[STIX]{x1D641}=\frac{\partial (x_{p},y_{p})}{\partial (x_{o},y_{o})}=\left(\begin{array}{@{}cc@{}}\displaystyle f^{\prime }\cos ^{2}{\it\theta}+\frac{r_{p}}{r_{o}}\sin ^{2}{\it\theta} & \displaystyle \left(f^{\prime }-\frac{r_{p}}{r_{o}}\right)\sin {\it\theta}\cos {\it\theta}\\ \displaystyle \left(f^{\prime }-\frac{r_{p}}{r_{o}}\right)\sin {\it\theta}\cos {\it\theta} & \displaystyle f^{\prime }\sin ^{2}{\it\theta}+\frac{r_{p}}{r_{o}}\cos ^{2}{\it\theta}\end{array}\right),\end{eqnarray}$$
in which
$|\unicode[STIX]{x1D641}|=r_{p}f^{\prime }/r_{o}$
. Therefore, the water depth tensor and gravitational acceleration in Cartesian coordinates become
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\unicode[STIX]{x1D65D}}=\frac{\unicode[STIX]{x1D641}\unicode[STIX]{x1D641}^{\text{T}}}{|\unicode[STIX]{x1D641}|}h=\frac{r_{o}}{r_{p}f^{\prime }}\left(\begin{array}{@{}cc@{}}\displaystyle f^{\prime 2}\cos ^{2}{\it\theta}+\frac{r_{p}^{2}}{r_{o}^{2}}\sin ^{2}{\it\theta} & \displaystyle \left(f^{\prime 2}-\frac{r_{p}^{2}}{r_{o}^{2}}\right)\sin {\it\theta}\cos {\it\theta}\\ \displaystyle \left(f^{\prime 2}-\frac{r_{p}^{2}}{r_{o}^{2}}\right)\sin {\it\theta}\cos {\it\theta} & \displaystyle f^{\prime 2}\sin ^{2}{\it\theta}+\frac{r_{p}^{2}}{r_{o}^{2}}\cos ^{2}{\it\theta}\end{array}\right)h, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{g}=|\unicode[STIX]{x1D641}|g=\frac{f^{\prime }r_{p}}{r_{o}}g, & \displaystyle\end{eqnarray}$$
$\tilde{\unicode[STIX]{x1D65D}}$
has eigenvectors of
$f^{\prime }r_{o}/r_{p}$
and
$r_{p}/(f^{\prime }r_{o})$
in the radial and azimuthal directions. It is easy to see that (A 4a
) and (A 4b
) are the same as (4.4), which was obtained using (2.10). The condition for a perfect physically realizable cloak (non-magnetic in the context of electromagnetic cloaking) is therefore obtained as
$|\unicode[STIX]{x1D641}|=C$
.
