Numerical and experimental study of an invisibility carpet in a water channel

We propose a numerical and an experimental study of an invisibility carpet for linear water waves. In the first part, we introduce the concept of an invisibility carpet in the case of linear water waves and apply this concept for a bounded problem: a wavetank. In the second part, we study a simpler case where we attempt to render invisible a vertical dihedral at the end of a wavetank. This is done by placing a structure consisting of 18 vertical poles with trapezoidal cross-sections in front of the dihedral. For these two configurations, with and without the carpet, we focus on the far-field reflected wave consisting of an inline mode and the first sloshing (plus progressive) mode. We show that our design achieves invisibility.


I. INTRODUCTION
In 2006, after some theoretical proposal on electromagnetic transparency with a plasmonic coating [1], the physicists Pendry, Schurig, and Smith published an original work theorizing that a finite-size object surrounded by a spherical coating consisting of a metamaterial defined by geometric transform might become invisible for electromagnetic waves [2].At the same time, Leonhardt proposed an analogous route to invisibility using conformal mappings [3].These three proposals have fuelled a new field of investigations in the area of metamaterials: "cloaking."The underlying idea behind cloaking using transformation optics is to map a point in optical space onto a spherical (invisibility) region, but it appears to be a severe limitation in the design of invisibility cloaks via transformation optics due to the singular behavior of the material parameters at the cloaks inner boundary.An alternative route is to use a one-to-one mapping to design an invisibility carpet, which is the bottom line of the bold proposal by Li and Pendry to conceal an object placed under a curved reflecting surface by imitating the reflection of a flat surface [4].Since the introduction of these concepts, we have seen some research work completing and extending the first ideas [5][6][7][8], with experimental validations [9][10][11][12][13].From electromagnetics, the topic has been extended to other fields such as acoustics [14,15], structural mechanics [16][17][18][19], plasmonics [20][21][22][23][24][25], and hydrodynamics [26][27][28][29][30].
In the first part of the paper, we develop concepts behind invisibility carpet with some numerical illustrations.In the second part of the paper we present a study, where a vertical dihedral, at the end of a wavetank, is made nearly invisible, in the framework of linearized potential flow theory.This study is rooted in the preliminary work presented in Ref. [31].

A. Basics
We first note it is possible to construct two linear operators defined on two different domains, yet sharing the same spectrum using a geometric transform.However, one of these * guillaume.dupont@fresnel.frtwo operators will necessarily have spatially varying, matrix valued coefficients.This mathematical property can be used in the design of metamaterials, whereby two different domains with Neumann boundary conditions behave in the same electromagnetic, acoustic, or hydrodynamic way (mimetism).To illustrate this property, we give an example of two distinct polygonal cavities sharing the same spectrum, and we further design an invisibility carpet for linear surface liquid waves in a curved channel.This structured metamaterial bends surface waves over a finite interval of Hertz frequencies.
We consider the following spectral problem: In the case of Neumann boundary conditions, the resolvent of operator A 1 is compact and its spectrum consists only of a countable set of discrete eigenvalues with a single accumulation point (0 or infinity depending upon whether we look at the operator or its inverse) [32].The question is: can we construct another operator A 2 , acting on functions in a domain 2 ( 1 and 2 are different domains), where the spectrum σ (A 2 ) is identical to σ (A 1 ).And the answer is positive, we can construct A 2 if we use a change of variables mapping the domain 1 on the domain 2 , cf.Fig. 1.The eigenvalues β and associated eigenfunctions ψ are where the eigenvalues β are real and positive (T is symmetric) and can be related (one by one) with the eigenvalues λ and eigenfunctions φ of Eq. (1).We numerically illustrate this fact with finite-element computations in Fig. 2, where two types of polygonal cavities mapped onto one another share the same spectrum.

B. Setup of the physical problem
Let denote the free surface of a channel occupied by a fluid, i.e., the interface between air and fluid; see Fig. 3.If we assume that the fluid is incompressible and irrotational, we know that the velocity field derives from a potential, which under the hypothesis of small perturbations of the free interface separating the fluid with ambient atmosphere leads to the T is the upper left block of the inverse of the symmetric matrix T = J T J/det(J).Our proposal of cloaking is to asymptotically approach the spectrum σ (A 2 ) with a sequence of spectra σ (A η ) associated with structured media of typical heterogeneity size η, when η goes to zero.
Helmholtz equation [32]: with k the spectral parameter related to the angular wave frequency ω (measured in radians per unit second) via the can be attributed to the fact that 2 is filled with a heterogeneous anisotropic medium.One should note that the eigenfields (a)-(e) and (f)-(j) associated with nonzero eigenvalues (arranged by increasing values), share the same features (such as mirror symmetry with respect to the y axis), up to a stretch along the y axis.Here, 1 is a triangle with vertices (−1,0), (1,0), and (2,0), and 2 is a polygone with vertices (−1,0), (0,1), (1,0), and (2,0).We map 1 onto 2 using Eq. ( 5) with y 1 (x) = a 1 x + b 1 and y 2 (x) = a 2 x + b 2 , where 6), we deduce the coefficients of the transformation matrix T: 4 }α, with α = (y 2 − y 1 )/y 2 .Importantly, T is a symmetric positive definite matrix; see Eq. ( 7).dispersion relation: Here, h is the depth of water in the channel, σ the surface tension at the free surface, ρ the fluid's density, and g the gravity.
Importantly, Eq. ( 3) also models the propagation of pressure waves in a compressible fluid, in which case k = ω √ ρ/λ, where ω is the angular pressure wave frequency, ρ is the unperturbed fluid mass density (a mass in kilogrammes per unit volume in metres cubed), and λ is the fluid bulk modulus (i.e., it measures the substances resistance to uniform compression and is defined as the pressure increase needed to cause a given relative decrease in volume, with physical units in pascals).
The linearized problem, Eqs. ( 3) and ( 4), allows for straightforward analogies between transverse electromagnetic and acoustic waves propagating in structured cylindrical domains, see Ref. [26] for the design of an invisibility cloak for surface liquid waves, theoretically and experimentally shown to work between 10 and 15 Hz for a structured cloak of radius 10 cm immersed in a fluid of depth 0.9 cm, with small surface tension and large density (σ = 13.6 N/cm and ρ = 1.529 kg/L, i.e., methoxynonafluorobutane), so that ω ∼ k √ gh.Here, we consider a channel of side lengths 1500 and 100 cm filled with a depth of water of 10 cm (σ ∼ 70 N/cm, ρ ∼ 1 kg/L), hence working frequencies should be scaled down roughly by a factor 10 with respect to Ref. [26], keeping in mind the dispersion relation Eq. ( 4) should play a more prominent role.

C. Design of a heterogeneous anisotropic fluid
Let us first introduce a simple geometric transform mapping the first region onto the second one.The bottom line is the bold proposal by Li and Pendry to conceal an object that is placed under a curved reflecting surface by imitating the reflection of a flat surface [4] in the context of electromagnetic waves in open space.In the present hydrodynamic case, the domain is bounded (water channel), its walls are rigid (Neumann conditions), and the geometric transform reads as follows: where J xx is the associated Jacobian Matrix, and α = (y 2 − y 1 )/y 1 .The effect of the transform is a stretch along the y axis,  4); (a) Field inside a straight channel filled with a homogeneous isotropic fluid corresponding to the Laplace operator A 1 in Eq. ( 1); (b) Field inside a curved channel filled with a homogeneous isotropic fluid for the same operator A 1 (note however the domain is different); (c) Field inside a curved channel filled with a heterogeneous anisotropic fluid described by Eq. ( 6) corresponding to the perturbed Laplace operator A 2 in Eq. ( 2).The color scale is in arbitrary units.The strong similarity between fields in (a) and (c) is noted.
which is parallel to the direction of wave propagation in the channel, i.e., the horizontal direction in Figs. 3 and 4.
The metric tensor associated with the transformed coordinates takes the following form (and its effect on the Cartesian metric is shown in Fig. 5, left panel, where the y axis is parallel to the vertical direction): It has been shown in Ref. [26] that the Block diagonal part of T −1 can be interpreted in terms of an effective shear viscosity associated with an anisotropic fluid.Indeed, the Helmholtz equation, Eq. ( 3), governing the propagation of linear surface water waves at the free interface, should be recast as Eq.
(2) after the geometric transform.However, the dispersion relation, Eq. ( 4), remains valid.We adopt the same viewpoint in the present note.It is interesting to look at the expression of the eigenvalues of T −1 as these are the relevant quantities to FIG. 5. Metrics associated with the Cartesian coordinate system (original domain, leftmost panel) and the transformed coordinate system (invisibility carpet, right panel) mapped onto one another via the transformation matrix T (note that the right angles are not preserved, i.e., the transformation is not conformal).design a structured channel: We note that λ 1 and λ i , i = 2,3, are strictly positive functions as obviously 1 2 and also α > 0. This establishes that T −1 is not a singular matrix for a two-dimensional carpet, which is a big advantage over two-dimensional cloaks obtained by blowing up a point onto a disk [4]: The transformation matrix is indeed singular at an invisibility cloak's inner boundary (one eigenvalue goes to infinity, while the other two go to zero [33]).The physical interpretation for the lack of singularity in invisibility carpets' parameters is rather intuitive: in order to flatten the wavefront of a wave incident upon a bump, one has to decelerate the wave inside the carpet, and this can be easily achieved with conventional materials.On the contrary, in order to flatten the wavefront of a wave incident upon an invisibility cloak, one has to accelerate the wave inside the cloak, in such a way that there is no phase shift observed in the foreward scattering.In electromagnetism, this means that one requires the cloak to be superluminal, which leads to a time paradox (causality problem).This is in essence what the singularity in invisibility cloaks' parameters leads to.From a numerical point of view, it is also important to end up with a nonsingular T matrix in the case of invisibility carpets: The variational form associated with Eq. ( 2) has the classical properties for existence and uniqueness of the solution (symmetry, coercivity, and boundedness of the bilinear form), hence the finiteelement formulation can be straightforwardly implemented.In this article, we used the commercial package COMSOL MULTIPHYSICS.

D. Structured fluid
Let us now mimic the heterogeneous anisotropic fluid described by the matrix T −1 (shear viscosity) in Eq. ( 6) using an effective medium approach whereby an assembly of rigid cylinders judiciously located is now fixed to the bottom of the channel, as shown in Fig. 3.Such a design approximates the transformed metric shown in Fig. 5 and it creates both the required anisotropy (driven by the aspect ratio of inclusions) and heterogeneity (all inclusions are not the same, the structured medium is quasiperiodic).It is clear that such a design will only work to certain extent and moreover will be constrained by the working eigenfrequency (the larger the eigenvalue of the operator, the larger the discrepancy between the ideal and approximated cases, as the cross-sectional size of rigid cylinders should be small compared to the wavelength in the homogenization setting).In Figs. 6 and 7, we show some representative fields corresponding to given eigenfrequencies in the range 1.72 Hz < ν < 2.33 Hz for a curved channel with a carpet (Fig. 6) and without a carpet (Fig. 7).We emphasize that the wavefront of the fields is nearly flat in Fig. 6: Cloaking is therefore evidenced.These numerical results clearly show the positive effect of the structured carpet.

III. EXPERIMENTAL SETUP
For the experimental part of the study, we use the wave flume at Ecole Centrale Marseille, which is L = 17 m long and b = 65 cm wide.Note that from now on we use coordinates (x,y) with x directed along the water channel, i.e., we rotate the coordinate axes of the preceding part through an angle π/2.The water depth is set at h = 40 cm.The configuration, which is the basis of the study, consists of a rigid vertical plate installed at the end of the tank, from wall to wall, at an angle of 60 degrees, thereby achieving a dihedral [Fig.8(a)].
The second configuration is an invisibility carpet, consisting of 18 vertical poles, with trapezoidal cross-sections, placed in  front of the dihedral [Fig.8(b)].For these two configurations, we run series of regular waves with wavenumber k in the range (0.7 π )/b to (2.3 π )/b according to the capabilities of the wavemaker.The reflected wave system, in the far-field, is the superposition of two modes: The inline mode and the first sloshing (plus progressive) mode.To access and separate these two modes, we use an array of 5 wave gauges over the width of the tank, set at different inline positions [Fig.8(a)], meaning that the same experimental case is run as many times as different positions are used.

A. Numerical study of the dihedral
In the case of the dihedral taken on its own, we use a semianalytical method, described below, to validate a finiteelement model formulated within the framework of linearized potential flow theory, and solved numerically with the COM-SOL Multiphysics software.In our wall-sided geometry, the linearized velocity potential writes where a i is the amplitude of the incident waves, h the water depth, g the gravitational acceleration, and the angular frequency ω is linked to the wavenumber k by the dispersion relation: The reduced potential ϕ(x,y) satisfies the Helmholtz equation in the fluid domain, no-flow conditions at the solid walls and appropriate ingoing and outgoing conditions at x → −∞ (Fig. 8).
The geometry at the end of the tank consists of two overlapping subdomains.The first domain is the angular sector 0 R 2d and 2π/3 θ π with d = b √ 3/3 (inside the dashed blue contour in Fig. 9) where the reduced potential ϕ takes the form with J 3m the Bessel function of the first kind.The second domain is the semiinfinite strip −∞ < x d and 0 y b (inside the solid red contour on Fig. 9), where the reduced potential can be written as with λ n = nπ/b, and N 2 the largest integer n such that λ n < k.Enforcing that the two expressions of the reduced potentials ϕ 1 and ϕ 2 coincide in the common region (green dots), we can determine the unknown coefficients A m and B n .To solve this problem, the series are truncated at orders M, N , with a number of control points N pt (N pt M + N + 2) distributed over the common region (green points in Fig. 9) and the following quantity, is minimized.We represent in Fig. 10 the results for the coefficients |B 0 | and |B 1 | (solid line).We can see that for kb/π 1 only the inline mode exists and for kb/π > 1 a first sloshing mode exists: kb/π = 1 corresponds to the cut-off frequency.We note that a new sloshing (plus propagative) mode appears for all integer values of kb/π.
The results of finite-element computations appear in dotted values in Fig. 10 and perfectly match the semianalytical values (solid curves).Thanks to this benchmark, we now dispose of a finite-element model to perform numerical investigations in the case of interactions of linear water waves with a vertical rigid wall in wave flume.We use this model below to perform our study of an invisibility carpet.

B. Invisibility carpet
The invisibility carpet is defined introducing a geometric transform in the same way as in Ref. [4] for electromagnetic waves, and constructed using an effective medium approach, whereby an assembly of rigid vertical poles, with trapezoidal cross-sections, is now fixed (with the water channel) in front of the dihedral.The carpet extends 1.5 m from the end point of the tank and the void fraction is close to 50%.In a first step, we want to determine the number of the vertical poles, keeping constant the extension and the void fraction of the carpet.The number of inclusions in the inline and transverse direction is varied in the finite-element computations, and the validity of the carpet is tested.To this end, an efficiency function is defined from the quantity where a number of reference abscissas x n are taken, from the edge of the carpet toward the wavemaker.F takes the value F 0 in the case with no carpet (dihedral alone) and the efficiency is defined as F /F 0 .An efficiency less than one means that the wave field tends to be unidirectional, and an efficiency equal to zero means that the dihedral has been made invisible.
The results of this optimization are presented in Fig. 11.We first note that the 18 inclusions configuration (yellow curve with square symbols in Fig. 11) appears to be the best, but it does not fully agree with effective medium theory, which predicts that the larger the number of poles and the smaller their cross-section, the lower E. The wall-sided geometry needs to be considered and could explain why a configuration with a low number of inclusions (18 inclusions) appears to work better than configurations with more inclusions (blue and red curves in Fig. 11).We choose the 18 inclusions' configuration for our invisibility carpet, and we numerically test this device before starting our experimental campaign.We report in Fig. 12 the comparison of the coefficients |B 0 | and |B 1 | of Eq. ( 12), from COMSOL computations, between the dihedral alone and the carpet placed in front of the dihedral.We note that the coefficient |B 0 | is close to 1 for the carpet [dashed red (dark gray) curve in Fig. 12], meaning that the inline mode dominates in the reflected wave system.However, the sloshing mode exists but is strongly decreased as compared to the dihedral alone [dashed blue (light gray) curve in Fig. 12].To illustrate the behavior of the reflected wave system, we present in Fig. 13 the calculated patterns of the free surface elevation for three nomalized wavenumber kb values with and without invisibility carpet.The first value (kb/π = 0.635) is below the cutoff frequency, hence the wave pattern is unidirectional in both cases.And, as shown in Fig. 12, for the other two values of kb/π above the cutoff frequency, the wave pattern is nearly unidirectional, but we can see a weak contribution of the sloshing mode.
Based on the numerical study, we build the invisibility carpet with 18 inclusions, and we place it at the end of the wavetank, in front of the dihedral.We run a series of regular waves (generated with a wavemaker) in the range (0.7 π )/b to (2.3 π )/b, with and without the carpet.Figure 14 shows the experimental |B 0 | and |B 1 | coefficients, as derived from the wave gauges measurements.The agreement with the results in Fig. 10 is excellent in the dihedral-alone case.One notes the slight oscillation for the B 1 coefficient in the dihedral and carpet cases, which can be attributed to genuine experimental inaccuracies.Importantly, when one  adds the carpet, the B 0 coefficient always takes higher values than for the dihedral-alone case, which is consistent with the fact that the B 1 coefficient is always lower: This demonstrates that the phase of the reflected wave is restored when we add the carpet.

IV. CONCLUSION
In this article, we have conclusively shown using numerical simulations and experiments that one can control the surface wave reflection with a 1-meter scale structured carpet.In a first step, we validated a finite-element subroutine making use of the COMSOL MULTIPHYSICS package (by comparison with semianalytical calculus presented in the paper) in order to perform a numerical study for interactions between linear water waves and an invisibility carpet.In a second step, we carried out an experimental campaign to demonstrate the successful functionality of our carpet.We noticed the dissipative behavior of such a device.A further work consists in identifying and quantifying damping sources and in this way, a particle image velocimetry study in a parallel plane close to the free surface (in the fluid domain) has been conclusively performed.This large-scale acoustic metamaterial might pave the way for nonovertopping dykes.

FIG. 1 .
FIG. 1. (Color online) If 1 and 2 are two bounded domains that can be mapped onto one another, via a change of coordinates described by the Jacobian matrix J, two self-adjoint bounded operators A 1 = − and A 2 = −T −1 33 ∇ • T −1 T ∇, respectively, defined on L 2 ( 1 ) and L 2 ( 2 ) have identical spectra σ (A 1 ) = σ (A 2 ), where T −1T is the upper left block of the inverse of the symmetric matrix T = J T J/det(J).Our proposal of cloaking is to asymptotically approach the spectrum σ (A 2 ) with a sequence of spectra σ (A η ) associated with structured media of typical heterogeneity size η, when η goes to zero.

FIG. 3 .
FIG. 3. (Color online) Schematic diagram of the waterwave experiment: A plane incident wave propagates from left to right at the free surface between fluid and air.It is reflected by an array of rigid pillars immersed in the fluid (the invisibility carpet).

FIG. 4 .
FIG. 4. (Color online) Numerical simulations for a surface water wave at frequency ν = ω/2π = 1.99 Hz, which is linked to the eigenvalue k of the Laplace operator in Eq. (3) via Eq.(4); (a) Field inside a straight channel filled with a homogeneous isotropic fluid corresponding to the Laplace operator A 1 in Eq. (1); (b) Field inside a curved channel filled with a homogeneous isotropic fluid for the same operator A 1 (note however the domain is different); (c) Field inside a curved channel filled with a heterogeneous anisotropic fluid described by Eq. (6) corresponding to the perturbed Laplace operator A 2 in Eq. (2).The color scale is in arbitrary units.The strong similarity between fields in (a) and (c) is noted.

FIG. 6 .
FIG. 6. (Color online) Eigenfields for 1.72 Hz < ν < 2.23 Hz in a curved channel with the structured carpet.The flat wavefronts of all eigenfields is noted.

FIG. 7 .
FIG. 7. (Color online) Eigenfields for 1.72 Hz < ν < 2.23 Hz in the same curved channel as in Fig. 6 but without the structured carpet.The disturbed wavefronts for the eigenfields is noted (except for the first two leftmost eigenfields).

FIG. 8 .
FIG. 8. (Color online) Schematic representation of the two studied configurations where the waves propagate from the left to the right.(a) Dihedral configuration.(b) Carpet configuration.(c) Photography of our experimental carpet.

FIG. 9 .
FIG. 9. (Color online) Geometry at the end of the tank.Subdomains 1 (highlighted by dashed blue lines) and 2 (highlighted by solid red lines) and control points (green dots).

FIG. 11 .
FIG. 11. (Color online) Efficiency function F /F 0 for 3 numbers of inclusion and 4 different configurations.Best results are achieved for the carpet with 18 inclusions.

FIG. 13 .
FIG. 13. (Color online) Calculated wave patterns in the end part of the tank for 3 normalized wavenumber kb values, with (left panels) and without (right panels) the carpet.

FIG. 14
FIG. 14. (Color online) Results from experiments: Coefficients |B 0 | (green and brown curves or dark gray curves) and |B 1 | (blue and yellow curves or light gray curves) of Eq. (12) for configurations without (solid lines) and with (dashed lines) the carpet.