Controlling solid elastic waves with spherical cloaks

We propose a cloak for coupled shear and pressure waves in solids. Its elastic properties are deduced from a geometric transform that retains the form of Navier equations. The spherical shell is made of an anisotropic and heterogeneous medium described by an elasticity tensor C' (without the minor symmetries) which has 21 non-zero spatially varying coefficients in spherical coordinates. Although some entries of C, e.g. some with a radial subscript, and the density (a scalar radial function) vanish on the inner boundary of the cloak, this metamaterial exhibits less singularities than its cylindrical counterpart studied in [M. Brun, S. Guenneau, A.B. Movchan, Appl. Phys. Lett. 94, 061903 (2009).] In the latter work, C' suffered some infinite entries, unlike in our case. Finite element computations confirm that elastic waves are smoothly detoured around a spherical void without reflection.


I. INTRODUCTION
In 2006, Pendry, Schurig and Smith [1], and Leonhardt [2] independently proposed some design of invisibility cloaks for light. The same year, Milton, Briane and Willis laid the foundations of transformational physics [3]. All these works generated a lot of interest in the metamaterials community, notably in the design of electromagnetic [4], water wave [5] and acoustic [6][7][8][9][10] cloaks. The latter fuelled the interest in acoustic metamaterials [11], which enable a markedly enhanced control of pressure waves, including lensing, via artificial anisotropy [12].
In 2009, Brun et al. in [13], discussed the design of a cylindrical cloak for in-plane elastic waves with an asymmetric elasticity tensor, which was a modified version of the Willis' type [15] transformed equations derived in [3]. In 2011, these two works were encompassed by Norris and Shuvalov in a more general elasticity framework [16].
In parallel to the developments of metamaterials for bulk elastic waves, some theoretical [17,18] and experimental [19] progress was made in the control of flexural elastic waves in thin plates. In the case of thin plates, the transformed governing equations (e.g. Kirchhoff) have a simpler structure which helps engineer structured cloaks.
In the present letter, we investigate spherical cloaks for solid elastic waves using a radially symmetric linear geometric transform which depends upon a parameter. We discuss their underlying mechanism and illustrate the theory using a finite element approach which is adequate to solve the Navier equations in transformed anisotropic heterogeneous media with asymmetric elasticity tensors.

A. The equations of motion
The propagation of elastic waves is governed by the Navier equations. Assuming time harmonic exp(−iωt) dependence, with ω as the angular wave frequency and t the time variable, allows us to work directly in the spectral domain. Such dependence is assumed henceforth and suppressed, leading to where ρ is the density and σ the stress tensor of the (possibly heterogeneous isotropic) elastic medium, u = (u r , u θ , u φ ) and C are respectively the three-component displacement field and the rank-four (symmetric) elasticity tensor expressed in a spherical coordinate basis x = (r, θ, φ).

B. The transformed equations of motion
Let us consider the coordinate change x −→ x , where x = (r , θ , φ ) are stretched spherical coordinates. This leads to a transformed equation [16] arXiv:1403.1847v2 [cond-mat.mtrl-sci] 23 Apr 2014 in general. Here S and D are third order tensors possibly encompassing a iω dependence, ∇ is the gradient in transformed coordinates x and u (x ) = u (r , θ , φ ) is a transformed displacement in stretched spherical coordinates. Note that the transformed stress σ is generally not symmetric. Note also that in general A is a matrix field. One possibility to preserve the symmetry of the stress tensor, is to assume that A is a multiple A = ξ∂x /∂x, of the Jacobian matrix ∂x /∂x of the transformation, where ξ is a non-zero scalar, in which case one obtains a Willis-type equation [3,15]. However, there is one special case for which u = u when A is the identity matrix I, which leads to [13] where the body force is assumed to be zero, the elasticity tensor C does not have the minor symmetries (Cosserat Material) and the stretched density ρ is a scalar field. This equation is derived from (1) by noting that S = 0 when A = I, [16]. In the sequel, we work in the framework of (3). Let us consider the geometric transform which maps a sphere of radius 0 < r ≤ r 2 onto on shell r 1 < r ≤ r 2 , see Figure 1.
Design of transformation-based Cosserat elastic cloaks has been first discussed in the cylindrical case in [13] where it only involved a tensor C with 8 non-vanishing coefficients, whereas in the present spherical case, we need to consider a tensor C with 21 non-vanishing coefficients. Moreover, the displacement field has three components in our case.
, where r1 and r2 are the inner and outer radii of the spherical cloak, respectively. A sphere (left) with an isotropic homogeneous symmetric elastic constitutive tensor C and homogeneous scalar density ρ, is mapped onto a hollow sphere (right) with a heterogeneous asymmetric elastic constitutive tensor C and a heterogeneous scalar density ρ .
By application of transformation (4), in the region 0 < r ≤ r 2 , the Navier equations (1) are mapped into the equations (3) with where a = r2 r2−r1 and b(r ) = (r −r1) r a, and the elasticity tensor C has 21 non-zero spherical components, namely C r r θθ = C θθr r = C r r φφ = C φφr r = λb(r ), C. Singularity at the cloak's inner boundary As discussed above, one should note that the minor symmetries are broken and b(r 1 ) = 0 which means C θθθθ /C r r r r = C φφφφ /C r r r r is infinite at the inner boundary of the cloak (infinite anisotropy). Moreover, the off-diagonal components are constant or vanish at the boundary r = r 1 .
Physically, this means that shear and pressure waves propagate much faster in the azimuthal and elevation directions than in the radial direction on the surface of the inner boundary. This should result in a vanishing phase shift between an elastic wave propagating in an isotropic homogeneous elastic medium, and another one propagating around the concealed region.
One should also note that there are no infinite entries within the elasticity tensor. This is a feature of the 3D elastodynamic cloak which is less singular than its 2D counterpart (the cylindrical elastic cloak in [13] had a vanishing density on the inner boundary of the cloak as well, but some entries of the elasticity tensor going to infinity at the inner boundary).

D. Impedance matching at the cloak's outer boundary
Let us now note that when r = r 2 , the geometric transform (4) leads to r = r 2 . In this case, the transformed density is ρ = aρ as b(r 2 ) = 1 and the diagonal components of the transformed elasticity tensor reduce to C r r r r = (λ + 2µ)/a , C θθθθ = C φφφφ = (λ + 2µ)a , hence C r r r r = C rrrr /a and C θθθθ = aC θθθθ and C φφφφ = aC φφφφ on the outer boundary of the cloak. This expresses the fact that a stretch along the radial direction is compensated by a contraction along the azimuthal and elevation directions in such a way that elastic media (cloak and surrounding isotropic elastic medium) are impedance-matched at r = r 2 . The cloak's outer boundary therefore behaves in many ways as an impedance matched 'thin' elastic layer. However, we note that the components of C pose no limitations on the applied frequency ω from low to high frequency, unlike for the case of coated cylinders studied back in 1998 in the context of elastic neutrality by Bigoni et al. [20].
The fact that C does not depend on ω, i.e. the cloak consists of a non-dispersive elastic medium, makes it work at all frequencies, but one should keep in mind that any structured medium designed to approximate the ideal cloak's parameters (e.g. via homogenization) would necessarily involve some dispersion, and thus limit the interval of frequencies over which the cloak can work. Such a feature has been already observed in [18,19], for cloaking of flexural waves in thin-elastic plates.

III. NUMERICAL ILLUSTRATION
We would like now to numerically test the cloaking efficiency. For this, we implement the 3 4 spatially varying entries of the transformed tensor in Cartesian coordinates in the finite element package COMSOL MULTIPHYSICS. We mesh the computational domain using 1105932 tetrahedral elements, 36492 triangular elements, 1128 edge elements and 25 vertex elements, see Figure 2. This domain consists of an isotropic homogeneous elastic medium within a sphere of radius r 3 = 10 m, containing a void surrounded by the cloak (a heterogeneous, anisotropic elastic metamaterial in a spherical shell of inner radius r 1 = 2 m and outer radius r 2 = 4 m ) and a point force located at (4 m, −7 m, 0 m), where (0, 0, 0) is the center of the cloak. The sphere of radius r 3 = 10 m is itself surrounded by a spherical shell of inner radius r 3 and outer radius r 4 = 12 m, which is filled with an anisotropic heterogeneous aborptive medium acting as a (reflectionless) perfectly matched layer (PML). This elastic PML is deduced from a geometric transform where we consider the radial function s r (r) = 1 − i. This transform leads, in the same way as (4) did for C and ρ , to a heterogeneous anisotropic asymmetric elasticity tensor C and a scalar (homogeneous) density ρ , whose expressions are obtained by taking a = 1 − i, b(r ) = 1 in (5)- (6). Since the elastic waves are damped inside the PML and reach the outer boundary of the shell with a vanishing amplitude, we set either clamped or traction free boundary conditions at r = r 4 (we have checked this does not affect the numerical result). One should keep in mind that any numerical implementation in a finite element package requires a Cartesian coordinate system. In our case, we used COMSOL where the transformed elastic tensors C for the cloak and C for the PML had up to 81 non-vanishing spatially varying entries. A good way to detect any flaw in the numerical implementation is to compare the gradient of the solution to the problem for a time harmonic point source in an isotropic homogeneous medium (supplied with spherical elastic PML), see middle panel in Figures 3-4, and the solution to the same problem when we have a cloak surrounding certain stress-free spherical region (e.g. a void in soil) which amounts to assuming that σ · n = (C : ∇ u) · n = 0 where n is the outward unit normal to the boundary of the void, see lower panel in Figures 3-4. One can see that the deformation of the elastic medium outside the cloak is nearly identical to that of the isotropic homogeneous elastic medium (we use the normalized density ρ = 1 and Lamé parameters µ = 1 and λ = 2.3, for respectively the shear modulus and the compressibility). By comparison, the deformation of the elastic medium is clearly visible for the void when it is not surrounded by the cloak, see upper panel in Figures 3-4. The small discrepancy between the middle and lower panels in Figures 3-4 is attributed to the numerical approximation of the infinitely anisotropic tensor at the inner boundary of the cloak (one way to avoid the singularity would be to use Kohn's transfrom [14]). Plots of elastic deformation can be somewhat misleading, and we therefore add plots of magnitude of the elastic displacement field in Figure 5, where it should be noticed that the shaded region behind the void in the upper panel (drop of wave amplitude and phase shift) , is almost completely removed in the lower panel thanks to the cloak. Upon inspection of the middle and lower panels, one can clearly see that the elastic field scattered by a void clothed with the cloak is virtually indistinguishable from bare isotropic elastic space.

IV. SOME POSSIBLE APPLICATIONS AND CONCLUSION
Finally, we would like to stress that another aspect of the spherical elastic cloak which we designed is its capability to protect any object placed inside the spherical void within the cloak from incoming elastic waves. Applications in anti-earthquake devices make this feature quite interesting to investigate, as typical frequencies of earthquakes are from 0.1 to 10 Hz, which is compatible with the frequency 3 Hz for a void of radius 2 m studied in the numerical illustrations of this letter. Besides from that, one can easily scale up the void e.g. a void of radius 20 m surrounded by a spherical cloak of outer radius 40 m would react in exactly the same way to a concentrated load of polarization (1, 1, 1) with pulsating frequency 0.3 Hz located at a distance r = 80.62 m from the center of the cloak.
Moreover, the material parameters of our cloak are not frequency dependent, and we numerically checked that similar results to those shown in the present letter occur for frequencies from 1 Hz to 10 Hz (we computed ten evenly spaced frequencies): we were limited below this range by the accuracy of the spherical PML (the larger the wave wavelength, the more absorption needed in PML) and above this range we are limited by the computational resources (numerical computation at 10 Hz required about 2 million tetrahedral elements for a converged result). The same remark holds for a void of radius r 1 = 2 m surrounded by a cloak of outer radius r 2 = 4 m as in this case the range of Figure 4: Deformation of an isotropic elastic medium (middle panel), a stress-free spherical obstacle of radius 2 m (upper panel) and the same void surrounded by an elastic cloak (lower panel) of inner radius r1 = 2 m and outer radius r2 = 4 m, subjected to a concentrated load of polarization (1, 1, 1) with frequency 3 Hz located at a distance 8.062 m from the origin. Deformation of components ε11 = ∂u 1 ∂x 1 (left column), ε22 = ∂u 2 ∂x 2 (middle column) and ε33 = ∂u 3 ∂x 3 (right column) of the strain tensor ε.
frequencies where the seismic signal would be detoured around the void is from 0.1 Hz to 1 Hz. Therefore, our study might find some applications in seismic metamaterials [21,22].
We finally note that recent advances in fabrication and characterization of elastic metamaterials [23] could foster experiments in an approximate 3D elastic cloak. Of course, the metamaterial would only be able in practice to display a strong (not infinite) anisotropy on the cloak's inner boundary and it would only work throughout a finite range of frequencies. Its properties could be derived for instance from an effective medium approach in a similar way to what was proposed [18] and experimentally validated [18] for elastic waves in thin plates. However, detouring solid waves (even partially) around a void is much more challenging than what was done for surface waves in plates. We hope the prospect of civil engineering applications [22] will fuel the research in seismic cloaks.