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


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


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


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


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]


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]


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:


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:


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:


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


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


Integrating by parts, we obtain


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


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


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


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

2.3 Perfecly Matched Layers for flexural waves

Figure 1: Schematic diagram showing the problem to be modelled: The plane elastic wave is incident from left to right in the presence of a diffracting obstacle of arbitrary shape. The different domains forming the PMLs regions are denoted respectively , and .

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


whose Jacobian is given by


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:


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


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


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


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




Figure 2: (a) Decay of the plane wave when it penetrates in the domain denoted by or in Figure 1 we represent here the amplitude of the plane wave in the direction or . (b) Propagation of the plane wave along the direction when it lies in the central domain (left of the red vertical line) and when it lies in the PMLs domain showing clearly that its amplitude vanishes rapidly in .

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


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


We have


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


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

(a)     (b)

Figure 3: Right: Real part of the displacement field distribution in the vicinity of the cloaked elliptic rigid clamped obstacle. The source is located at point and its wavelength is , which is of the same order as the inner and outer radii (major and minor) of the cloak: , and , . Left: Real part of the displacement field scattered by a rigid clamped obstacle of minor radius m and major radius m for an incoming cylindrical wave of the same frequency.

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


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


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

Figure 4: A plane bending wave is incident from left to right onto a coating consisting of two concentric ellipses of the same dimensions as in Figure 3. The Young modulus and density of the shell are deduced from the transformation given in Equation (27). We can see that the wave fronts are rotated with an angle radians.

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:


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


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.

(a)     (b)


Figure 5: (a) Real part of the displacement field when the source is located inside the coating at the point . The field seems to be emitted by a shifted source located at point (b); (c) Real part of the displacement along the line for a source located inside the coating at point as shown in Figure 5 (a) (solid line) and for a source located in the homogeneous plate at point (dashed-starred line) deduced from the geometrical transform given in (26).


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.


  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).
Comments 0
Request Comment
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
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
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 description