Finite elements modelling of scattering problems for flexural waves in thin plates: Application to elliptic invisibility cloaks, rotators and the mirage effect

# Finite elements modelling of scattering problems for flexural waves in thin plates: Application to elliptic invisibility cloaks, rotators and the mirage effect

## Abstract

We propose a finite elements algorithm to solve a fourth order partial differential equation governing the propagation of time-harmonic bending waves in thin elastic plates. Specially designed perfectly matched layers are implemented to deal with the infinite extent of the plates. These are deduced from a geometric transform in the biharmonic equation. To numerically illustrate the power of elastodynamic transformations, we analyse the elastic response of an elliptic invisibility cloak surrounding a clamped obstacle in the presence of a cylindrical excitation i.e. a concentrated point force. Elliptic cloaking for flexural waves involves a density and an orthotropic Young’s modulus which depend on the radial and azimuthal positions, as deduced from a coordinates transformation for circular cloaks in the spirit of Pendry et al. [Science 312, 1780 (2006)], but with a further stretch of a coordinate axis. We find that a wave radiated by a concentrated point force located a couple of wavelengths away from the cloak is almost unperturbed in magnitude and in phase. However, when the point force lies within the coating, it seems to radiate from a shifted location. Finally, we emphasize the versatility of transformation elastodynamics with the design of an elliptic cloak which rotates the polarization of a flexural wave within its core.

Finite Elements; High-Order Differential Operator; Perfectly Matched Layers; Transformation Elastodynamics; Cloaking

M. Farhat]Mohamed Farhat\corrauth, Sébastien Guenneau and Stefan Enoch

\ams

52B10, 65D18, 68U05, 68U07

## 1 Introduction

In 2006, Pendry et al. have shown that by surrounding a finite size object with a coating consisting of a metamaterial, we can render this system transparent to electromagnetic radiation [1]. The cornerstone of this work is the geometric transformation

 ⎧⎪ ⎪⎨⎪ ⎪⎩r′=R1+r(R2−R1)/R2,0

where , and are radially contracted cylindrical coordinates and is the Cartesian basis. This transformation maps the disk onto an annulus . In other words, if a source located in radiates in vacuum, the electromagnetic field cannot reach the disk and therefore this region is a shelter for any object. In layman terms, a point in the space is transformed into a disc leading to the creation of a hole in the physical space where waves cannot penetrate from outside but are guided around this area in a way that anything we put inside this disc would be invisible to exterior observers, in the context of electromagnetism [1, 2, 3, 4, 5], or neutral, in the context of elastodynamics [6, 7, 8, 10, 11, 9, 12].

From a geometric point of view, the change of coordinates means that in the annulus we should work in a stretched space with associated metric tensor

 T=Jtrr′Jrr′/det(Jrr′), (2)

with the Jacobian of the transformation from free to contracted cylindrical coordinates, its transpose and its determinant.

In terms of material parameters, the only thing to do in the annulus is to replace the material (homogeneous and isotropic) by an equivalent one that is inhomogeneous (its characteristics are no longer piecewise constant but merely depend on , , coordinates) and anisotropic ones (tensorial nature) whose properties are given by [14, 15, 13, 3]

 ε′––––=εT−1,andμ′––––=μT−1. (3)

Variants of this formula appear in many instances in the existing literature, such as in the book on the geometry of electromagnetism by Post, back in 1962 [16]. However, it seems that Pendry and Ward were the first ones to apply it to simplify the computational analysis of certain types of photonic materials in 1996 [17].

The diagonalised form of the tensors and was given in [1]

 εr=μr=r−R1r,εθ=μθ=rr−R1,ε3=μ3=(R2R2−R1)2r−R1r (4)

Mimicking the heterogeneous and anisotropic nature of the tensors of permittivity and permeability given by eq. (4) became possible only when Pendry suggested the use of newly discovered metamaterials which are composite structures manufactured for their exotic properties enabling one to control the electromagnetic field. A team lead by Pendry and Smith implemented this idea using a metamaterial consisting of concentric layers of Split Ring Resonators (SRR), which made a copper cylinder invisible to an incident plane wave at GHz as predicted by the numerical simulations [2]. Independently, Leonhardt studied conformal invisibility by solving the Schrödinger equation which is valid in the geometric optics limit [18].

A different path to invisibility was followed by McPhedran et al. [19] who proposed to cloak a countable set of line sources using anomalous resonance when it lies in the close neighbourhood of a cylindrical coating filled with a negative index material which is nothing but a cylindrical version of the celebrated perfect lens of Sir John Pendry [20].

An impedance matching route to invisibility was further developed by Alu and Engheta [21]. Although very promising, their proposal relies on a specific knowledge of the material properties of the object being concealed. Recently, these authors gave the experimental proof of their idea [22].

Other routes to invisibility are based on homogenization; A team led by Shalaev [4] has shown the possibility to make an object nearly invisible in TE polarization. A reduced set of material parameters was introduced to relax the constraint on the permeability, thus leading to an impedance mismatch with vacuum which was shown to preserve the cloak effectiveness to a good extent.

Farhat et al. [23] analyzed cloaking of transverse electric (TE) fields through homogenization of radially symmetric metallic structures. The cloak consists of concentric layers cut into a large number of small infinitely conducting sectors which is equivalent to a highly anisotropic permittivity. This structure was shown to work for different wavelengths provided they are ten times larger than the outermost sectors.

Soon after, Cummer and Schurig [6] analyzed the acoustic cloaking for pressure waves in a transversely anisotropic fluid by exploiting the analogy with TE electromagnetic waves. Torrent and Sanchez-Dehesa [7] subsequently investigated this cloaking for concentric layers of solid lattices behaving as anisotropic fluids in the homogenization limit. Using a similar approach, Farhat et al. [8] independently demonstrated cloaking of surface liquid waves using a micro-structured metallic cloak which was experimentally validated at Hz. Quite remarkably, Chen and Wu [9] and Cummer et al. [24] noticed that a acoustic cloaking for pressure waves in a fluid can be envisaged since the wave equation retains its form under geometric changes.

However, when one moves to the domain of elastodynamics, Milton, Briane and Willis [10] have shown that the governing equations are not invariant under coordinate transformations and consequently that if cloaking exists for such type of waves, it would be of a different nature than its acoustic and electromagnetic counterparts. Indeed, the equation of propagation of elastic waves in solids (Navier equation) with a time harmonic dependence can be written in the weak form:

 ∇⋅C:∇u+ρ0ω2u=0, (5)

where , is the three-component vector displacement field, is the scalar density of the elastic medium, is the rank elasticity tensor, is the wave angular frequency, and is the time.
Moreover, Milton, Briane and Willis noticed [10] that by introducing the geometric transform such that with , thereby enforcing the symmetry of the elasticity tensor, equation (7) takes the form:

 ∇′⋅(C′+S′):∇′u′+ρ––––′ω2u′=D′:∇′u′, (6)

which importantly preserves the symmetry of the new elasticity tensor . However, this transformed equation contains and which are two rank 3 (symmetric) tensors such that and is a rank tensor whose expressions can be found in [10].

Nevertheless, Brun, Guenneau and Movchan recently discovered that if one does not constrain the symmetry of the elasticity tensor, equation (7) takes the form for the case of fully coupled in-plane shear and pressure waves:

 ∇⋅C′:∇u+ρ′ω2u=0, (7)

where is a rank four tensor with only the main symmetries and a scalar quantity, both of them spatially varying, see [11].

Farhat, Guenneau, Enoch and Movchan further considered the case of biharmonic bending waves propagating in thin plates and have shown that by considering a radially dependent isotropic mass density and a radially dependent and orthotropic flexural rigidity, it was possible to extend cloaking by refraction to the case of out-of-plane elastic waves [25].

In this paper, we give a numerical analysis of the biharmonic problem based on the Finite Elements Method (FEM) which is an approximation technique for the partial differential equations (PDE) and require a variational formulation of the problem to be studied (expressing expressions in a weak form). Indeed, the COMSOL formulation of the fourth order biharmonic equation is given, supplied with appropriate boundary conditions (clamped or stress-free) and Perfectly Matched Layers (PMLs). In the last section, elliptic cloaking is numerically performed and the rotator and mirage effects established for these waves confirming that they behave in a way similar to electromagnetic and acoustic waves.

## 2 Finite Element modelling for the biharmonic equation

The equations for bending of plates are well known and can be found in many textbooks, such as those of Timoshenko or Graff [26]-[27]. The displacement is in the -direction. We choose to work in cylindrical coordinates with a time harmonic dependence i.e. , where is the time and is the angular frequency.

The wavelength is supposed to be large enough compared to the thickness of the plate and small compared to its in-plane dimensions . In this case we can adopt the hypothesis of the theory of Von-Karman.

To derive the equation of motion of the flexural out-of-plane waves, we can use a variational point of view [28].

### 2.1 Weak form associated with the fourth-order equation

We now consider the following fourth order elliptic equation

 a∇⋅(b∇(a∇⋅(b∇W)))=finΩ, (8)

where is a scalar, and is a matrix describing the material of the domain, and a source term. This equation is supplied with ’Dirichlet’ (or clamped) boundary conditions and on (we note that the second boundary condition is needed since we are in presence of a fourth-order partial differential equation).

By multiplying Equation (8) with a test function and integrating on the domain

 ∫Ωdτ[(∇⋅(b∇(a∇⋅(b∇W)))W′]=∫ΩdτfW′. (9)

Integrating by parts, we obtain

 −∫Ωdτ[(a∇⋅(b∇W))(a∇⋅(b∇W))]=∫ΩdτfW, (10)

where we have used the fact that and its normal derivative vanish on .

Importantly, we need to check that this problem is well posed. For this, we denote by (resp. ) the (class of) square integrable functions on , with values in (resp. ). We also introduce the Hilbert space

 H3(Ω)={v∈L2(Ω):∇v∈L2(Ω),Δv∈L2(Ω),∇Δv∈L2(Ω)} (11)

If we further assume that and on , the above space is denoted .

The existence and uniqueness of the solution is then straightforwardly ensured by the Lax-Milgram theorem applied for , and , provided that there exist , , , such that and for all .

### 2.2 Discretisation of the variational equation

To discretise this variational problem, the basic idea is to approximate with a linear combination of simple functions called form functions and spanning a functional space of finite dimension (i.e. a Galerkin space with a finite number of degrees of freedom). Then, this development is inserted into the weak form of the fourth-order partial differential equation thus generating a finite system of equations.

It is common to consider a triangular mesh with a reference cell the discussion we assume a square mesh. We then consider e.g. a basis of second order polynomials which spans a space of dimension 6. High-order polynomial basis can also be used. The function can be written as

 W=N∑iβiϕi(x,y) (12)

where are values of at the node points of the triangular mesh. The form functions satisfy the condition ( is the Kronecker symbol).

Finally using the development of (12), we obtain the system of n equations with n unknowns where we replaced by the set of form functions

 −N∑j=1(∫Ωdτ[(a∇⋅(b∇ϕj))(a∇⋅(b∇ϕi))])=∫Ωdτfϕi=0,i=1,...,N, (13)

and this can be written as a linear algebraic system which can be solved iteratively.

### 2.3 Perfecly Matched Layers for flexural waves

Perfectly Matched Layers (PMLs) are equivalent to a geometrical transformation. In Cartesian coordinates they can be expressed as

 xs(x)=∫x0dx′sx(x′)ys(y)=∫y0dy′sy(y′)andzs(z)=∫z0dz′sz(z′), (14)

whose Jacobian is given by

 Jxxs=∂(x,y,z)∂(xs,ys,zs)={% \bf Diag}(1sx,1sy,1sz) (15)

where Diag represents a diagonal matrix.

We consider now the transformation of the biharmonic equation under an arbitrary change of coordinates (say from to ). To do so, we have to consider the way the two equations in the following system transform:

 ⎧⎨⎩∇.(−ζ–––−1∇W)+λ−1V=0∇.(−ζ–––−1∇V)+λ−1β40W=0 (16)

where is an inhomogeneous anisotropic tensor and is an inhomogeneous coefficient of the material ( and ).

For example, the first equation of this system is similar to Helmholtz’s equation except that here we have two functions and instead of a unique function. Though we can write it in an integral form in the domain of validity let’s say

 ∫Ωdτ∇.(−ζ–––−1∇W)+∫Ωdτλ−1V=0 (17)

with the infinitesimal element of volume.
By multiplying this equation by a test function and integrating by parts, we obtain

 −∫Ωdτ∇ψ⋅(−ζ–––−1∇W)+∫Ωdτλ−1ψV=0 (18)

Now, by operating the previously described coordinate transform characterized by its jacobian that we will will note simply to lighten the notations, (18) transforms as

 −∫Ω′dτ′1det(J)Jt∇′ψ⋅(−ζ–––−1Jt∇′W)+∫Ω′dτ′1det(J)λ−1ψV=0 (19)

where we have used the generic way the gradient and the infinitesimal volume transform: and . Finally, using the useful expression of the scalar product where is the transpose of we can write this equation in the more adapted form

 ∫Ω′dτ′(∇′ψ)t{Jζ–––−1Jtdet(J)}∇′W+∫Ω′dτ′{λ−1det(J)}ψV=0 (20)

Hence, we can see from this last expression, that transforming the biharmonic equation from the initial system of coordinates to the new one is equivalent to replacing the initial parameters and by the deduced parameters and

 {ζ′––––−1=Jxsxζ–––−1Jtxsx/det(Jxsx)λ′−1=λ−1/det(Jxsx) (21)

We can obtain the same expressions by considering the second equation in (16).

Now, consider that we have an isotropic and homogeneous media, the use of PMLs introduce the complex matrix which represents the inverse of the metric tensor. It can be expressed in terms of the parameters of the Jacobian of the PML transformation , and

 T−1PML={\bf Diag}(syszsx,sxszsy,sxsysz) (22)

We have

 sz=sy=1,sx(x)={1% In the domain of study1−iσyωεin the region R1(% See Figure ???) (23)

For the regions , we have just to replace by and 1 by 2 in equation (23), and for the regions , we have

 sz=1,sy(y)=1−iσyωε% andsx(x)=1−iσxωε (24)

Figure 2 clearly confirms the mechanism of PMLs based on the transformation of biharmonic equation from the initial system () to () representing absorbing layers in such a way incident waves are perfectly transmitted (without any reflexions) independently from their frequency or incidence angles.
To our knowledge, this is the first time PMLs were applied to the fouth order biharmonic problem.
For exapmle, Figure 2 (b) shows the decay of an incident plane wave in the PML zone (rectangle of Figure 1): the amplitude of the wave can be considered zero after only a short propagation distance in the absorber medium.

## 3 Elliptic invisibility cloaks and the mirage effect for biharmonic problems

### 3.1 Elliptic invisibility cloak for flexural waves

This may be deduced from the circular case by scaling the Cartesian coordinates. The global transformation is a mapping of a holey elliptic domain (the inner and outer boundaries are concentric ellipses with the same eccentricity) on a simply connected elliptic domain bounded by the outer ellipse of the cloak. The detailed computation of the equivalent material properties is given by the following sequence of transformations. The starting point is the ellipse that bounds the exterior limit of our cloak with its principal axes chosen conveniently parallel to the coordinates axes. The first step is aimed at restoring the previous situation namely a circular cloak. For this purpose, the plane is scaled by a factor along the (arbitrary chosen) -axis so that the initial ellipse becomes a circle (in the scaled coordinates) defined by the transformation (and simply ) characterized by the Jacobian matrix . The next three transformations are then the ones used to build the circular cloak: transformation to cylindrical coordinates, radial contraction (the active part), and transformation to rectangular coordinates. It is important to note that it is the scaled variable that is involved in these various operations. The last step to be performed is an inverse scaling along -axis: to recuperate the initial elliptical shape of the cloak and an identity for the transformation of the outside of the cloak. At the end, a cloak is obtained whose inner and outer boundaries are ellipses with -axis of the hole, , -axis of the hole, -axis of the external boundary, , -axis of the external boundary. The total Jacobian of this sequence of transformations is

 Jxx′=JxxsJxsrJrr′Jr′x′sJx′sx′ (25)

The inverse of the matrix is given explicitly by [5]

 T−1=Diag(1,sy,1)R(θ′s)Diag(α,r′sr,1)R(−θs)Diag(1,1s2y,1)R(θs)Diag(α,r′sr,1)R(−θ′s)Diag(1,sy,1)rαr′s, (26)

with , and .

Note the angles and distances are computed in the scaled coordinate systems where the ellipses are mapped to circles.

Figure 3 (b) shows the calculated displacement field distribution when an elastic point source vibrating in the out-of-plane direction is placed near the optimized cloak. It is shown that in the near region of the cloak, the displacement field stays almost unperturbed, thus a so called perfect elliptic cloak is obtained. In order to see how much improvement has been achieved through the presence of the elliptic coat, snapshots of the displacement in the presence of a clamped obstacle ( and where describes the boundary of the ellipse) are shown in (a).

### 3.2 Transformation elastodynamics: Elliptic rotator and mirage effect

Based on the coordinate transformation methods [18], we propose here a two-dimensional transformation-media coat that can rotate fields but remains itself invisible. The Young modulus tensor and density of such a ”field rotator” for flexural waves can be deduced from appropriate coordinates transforms.

The difference with invisibility (considered in the previous section and where we map a point on a circle of radius ) is that here, we map a circle (of radius or ) onto itself. So let’s consider the following transformation:

 ⎧⎪⎨⎪⎩r′=rθ′=θ+α1r+β1z′=z (27)

with and . is the angle of rotation of the wave fronts.
The Jacobian of this transformation is given by

 Jrr′=∂(r,θ,z)∂(r′,θ′,z′)=⎛⎜⎝100−α110001⎞⎟⎠ (28)

This expression may be injected in Equation (25) instead of the invisibility Jacobian to obtain the matrix which gives us the physical parameters of the coat (Young’s modulus and density). In addition, inside the inner circle of radius , we have the additional transform (, and ) whose Jacobian is the identity matrix and doesn’t contribute thus to the materials specifications inside the cloak. Finally, the FEM computations using these expression show in Figure 4 the displacement filed in the vicinity of an ideal field rotator designed by transformation-media theory.

We consider also another application of coordinates transformation and which has been verified for others types of waves (acoustic and elastodynamics).

In Figure 5 the vertical displacement outside the cloak seems to originate from a location , which is slightly shifted with respect to the real position of the source (located inside the cloak itself). This effect is similar to the mirage effect already observed in electromagnetic cloaks [3]. The vertical displacement on the line is shown for both a source located inside the coating and a source shifted laterally at in a homogeneous plate (i.e., without any cloak and inclusion). The respective positions are given in the figure caption.

In conclusion, we have studied in this paper some properties of the biharmonic equation governing the propagation of flexural waves in thin elastic plates. we have shown the way to generalize some concepts like cloaking or the rotator and the mirage effect to the domain of elasticity in the special case of elliptic geometries. Despite the fact that general elastic equations are not invariant under geometrical changes, we have demonstrated that cloaking through coordinates changes is applicable to some special configurations (thin plates). Moreover Perfectly Matched Layers have been performed using the finite elements technique for the biharmonic problems.

## Acknowledgments

The authors would like to thank insightful discussions with Prof. A.B. Movchan and N.V. Movchan at Liverpool University, Dr. C.G. Poulton at the University of Technology of Sydney (UTS) and Prof. R.C. McPhedran at the University of Sydney.

### References

1. J.B. Pendry, D. Schurig, and D.R. Smith, Science 312, 1780 (2006).
2. D. Schurig, J.J. Mock, B.J. Justice, S.A. Cummer, J.B. Pendry, A.F. Starr, D.R. Smith, Science 314, 977 (2006).
3. F. Zolla, S. Guenneau, A. Nicolet, and J.B. Pendry, Opt. Lett. 32, 1069-1071 (2007).
4. W. Cai, U.K. Chettiar, A.V. Kildiev and V.M. Shalaev, Nature 1, 224-227 (2007).
5. A. Nicolet, F. Zolla and S. Guenneau, IEEE Trans. MAg. 44, 1150-1153 (2008)
6. S.A. Cummer and D. Schurig, New J. Phys. 9, 45 (2007).
7. D. Torrent and J. Sanchez-Dehesa, New J. Phys. 10, 063015 (2008).
8. M. Farhat, S. Enoch, S. Guenneau and A.B. Movchan, Phys. Rev. Lett. 101, 134501 (2008).
9. H. Chen and C. T. Chan, Appl. Phys. Lett. 91, 183518 (2007).
10. G.W. Milton, M. Briane, and J.R. Willis, New J. Phys. 8, 248 (2006).
11. M. Brun, S. Guenneau and A.B. Movchan, Appl. Phys. Lett. 94 061903 (2009).
12. H. Chen, B.I. Wu, B. Zhang, and J.A. Kong, Phys. Rev. Lett. 99, 063903 (2007).
13. U. Leonhardt and T. G. Philbin, New J. Phys. 8, 247 (2006).
14. A. Nicolet, J.F. Remacle, B. Meys, A. Genon and W. Legros, Appl. Phys. 75, 6036-6038 (1994).
15. F. Zolla, G. Renversez, A. Nicolet, B. Kuhlmey, S. Guenneau and D. Felbacq, Foundations of photonic crystal fibres (Imperial College Press, London, 2005).
16. E.G. Post, Formal Structure of Electromagnetics; General Covariance and Electromagnetics (Interscience, 1962).
17. A.J. Ward and J.B. Pendry, J. Mod. Opt. 43, 773-793 (1996).
18. U. Leonhardt, Science 312 1777-1780 (2006).
19. N.A. Nicorovici, R.C. McPhedran and G.W. Milton, Phys. Rev. B 49, 8479-8482 (1994).
20. J.B. Pendry, Phys. Rev. Lett. 86, 3966-3969 (2000).
21. A. Alu and N. Engheta, Phys. Rev. E 95 016623 (2005).
22. B. Edwards, A. Alu, M. G. Silveirinha, and N. Engheta, Phys. Rev. Lett. 103 153901 (2009).
23. M. Farhat, S. Guenneau, A.B. Movchan and S. Enoch, Opt. Express 16, 5656-5661 (2008)
24. S.A. Cummer, B.I. Popa, D. Schurig, D.R. Smith, J. Pendry, M. Rahm, and A. Starr, Phys. Rev. Lett. 100, 024301 (2008).
25. M. Farhat, S. Enoch, S. Guenneau and A.B. Movchan, Phys. Rev. B 79 033102 (2009).
26. S. Timoshenko, Theory of plates and shells (McGraw-Hill, New York, 1940).
27. K.F. Graff, Wave motion in elastic solids (Dover, New York, 1975).
28. L.D. Landau, and E.M. Lifschitz, Elasticity theory (Pergamon Press, 1954).
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters

220662

How to quickly get a good answer:
• Keep your question short and to the point
• Check for grammar or spelling errors.
• Phrase it like a question
Test
Test description