Transformation cloaking and radial approximations for flexural waves in elastic plates

It is known that design of elastic cloaks is much more challenging than the design idea for acoustic cloaks, cloaks of electromagnetic waves or scalar problems of anti-plane shear. In this paper, we address fully the fourth-order problem and develop a model of a broadband invisibility cloak for channelling flexural waves in thin plates around finite inclusions. We also discuss an option to employ efficiently an elastic pre-stress and body forces to achieve such a result. An asymptotic derivation provides a rigorous link between the model in question and elastic wave propagation in thin solids. This is discussed in detail to show connection with non-symmetric formulations in vector elasticity studied in earlier work.


Introduction
There is a theoretical and practical interest in wave cloaking in the context of metamaterials, as outlined in [1][2][3][4][5][6][7][8]. Dynamic effects include anisotropy and localization [9,10], which can be be interpreted in the context of homogenization theory. In this regard, we would like to refer to the work [11][12][13] addressing the notion of an effective dynamic mass density in structured composites and acoustic materials, as well as analytical studies of dynamic localization in phononic crystals. The approach of dynamic homogenization has been systematically applied in [14,15] to vibrations of inertial lattice systems. The idea of a so-called 'geometric optics' transformation leading to a radially symmetric 'push-out' cloak is commonly used for computational and experimental implementation [1,[16][17][18][19]. In scalar problems, where the governing equation is reduced to Helmholtz form, such a transformation proves to be extremely efficient, leading to a model of a specially designed highly anisotropic inhomogeneous material occupying the cylindrical cloaking layer and channeling incident waves around a finite scatterer (an inclusion or a void). The continuum model of an invisibility cloak leads to singular behaviour of the theoretical material at the inner boundary of the cloaking region adjacent to the scatterer. In a practical implementation, a continuum cloak is replaced by a micro-structured composite, and examples of such implementation include water waves [8], flexural plates [18] and acoustics [20]. This micro-structure makes the cloak approximate, and such an approximation is frequency sensitive. A special challenge is presented for vector problems of elasticity discussed in [4,7,16,21].
The present paper addresses cloaking for flexural waves in Kirchhoff elastic plates. Firstly, we show that the governing equations are not invariant with respect to the radial 'push-out' transformation [3,22]. This observation implies that the cloaking design procedure, well developed for acoustics, vibration of elastic membranes and anti-plane shear problems (see, for example, [23,24]), does not apply to problems of flexural vibrations of elastic plates. Elastic Kirchhoff plates possess flexural rigidity and their out-of-plane vibrations are governed by a fourth-order partial differential equation. One of the main challenges appears to be the presence of propagating and evanescent waves representing solutions of the Helmholtz and modified Helmholtz equations, and the coupling of such waves via the boundary and interface contact conditions. In numerical simulations, it is apparent that in many configurations the flexural waves are led by their Helmholtz component (see, for example, [25,26]). However, for cloaking problems the multi-scale nature of a metamaterial makes the problem more challenging and it is not apparent that such decoupling is possible.
There is strong experimental evidence, as published in [18] and also outlined in [27], that within a predefined frequency range a by-pass system can be implemented around a finite obstacle in a flexural Kirchhoff plate. Such a by-pass system is evidently an approximate cloak, which would benefit strongly from an in-depth analysis paving the way to a broadening of the frequency range for the cloaking effect.
We explain the derivation of such an approximate cloaking model and present illustrative numerical examples which agree with the experimental evidence [18].
It was shown in [19], for a model of a square cloak, that a formulation for flexural waves in a Kirchhoff plate, after the cloaking transformation, includes additional terms in the governing equation; these may represent in-plane body forces and pre-stress. This approach provides a consistent procedure justifying the additional terms in the governing equation and cloaking is effective across the whole frequency range admissible for the plate model. Motivated by the results of [19], we also develop the full cloaking model for the radial 'push-out' transformation, and obtain an explicit closed-form representation for the pre-stress required to have a broadband cloak for flexural waves.
Finally, we present a detailed asymptotic analysis, which establishes a connection between the transformed equations for the fourth-order model of flexural waves and those for a vector problem of elasticity in thin solids.
2. Application of the radial 'push-out' transformation to a Kirchhoff-Love plate We begin with a simple case of the equation governing the out-of-plane displacement amplitude w X ( ) of an orthotropic homogeneous plate, in the absence of applied in-plane forces, under pure bending. As in [28], the fourth-order partial differential equation is where D R , Θ D and Θ D R are the flexural rigidities, ρ and h are the mass density per unit volume and thickness of the plate, respectively, and ω is the angular frequency.
The constitutive relations that define the moments are where ν R and ν Θ are the values of the Poisson ratios in the radial and tangential directions, respectively. We also note that , and the rigidities D R and Θ D satisfy the following symmetry relation:

R R
Further, if the plate is isotropic and homogeneous, then equation (1) is the flexural rigidity of the isotropic plate, so that the equation of motion (1) simplifies to The Jacobi matrix F in cylindrical coordinates θ r ( , ) has the form     4 for the normalized mass density, then it is tempting to assume that (9) represents an orthotropic inhomogenous plate with the stiffness rigidities (8). The question is: can such an assumption be justified?
On one hand, the fourth-order terms in (9) agree with the structure of (1). On the other hand, the additional lower-order terms have to be analysed.
We also note that, after the normalization, equation (7) can be written in the compact form where the differential operator  R We would like to emphasize that for > R 0 1 , the operators  R 2 1 and ∇ 2 are not the same. The operator  R 2 1 will be referred to as the 'shifted Laplace operator', which becomes the classical Laplace operator only when , i.e. in the absence of the cloak. Correspondingly, we will use the terms 'shifted Helmholtz' and 'shifted modified Helmholtz' for the operators β , respectively. Following the representation (10) we can express the transformed equation in the form  and hence w is the superposition of waves of shifted Helmholtz type w HS and shifted modified Helmholtz type w MS . A semi-analytical solution can be found by implementing the series representation In equation (17) J n is the Bessel function, H n (1) is the Hankel function, and I n and K n are the modified Bessel functions related to J n and H n (1) by respectively (see [29], equations 9.6.3 and 9.6.4). The coefficients of the expansion (16), (17) are determined from the boundary and the interface conditions on the contour of the cloak.

Transformation cloaking for a membrane versus flexural plate
The radial 'push-out' transformation (5) can be used to design a cloak that will route an incident wave around a finite-size obstacle in an elastic membrane. Norris [23] has discussed this problem in detail. The governing equation for the time-harmonic out-of-plane displacement u of an elastic membrane has the form where μ stands for the stiffness matrix and ρ and ω are the mass density and the radian frequency, respectively. If an invertible mapping  = x X ( ) is applied within the cloaking region, then the transformed equation becomes ⎛ It is important to note that equation (20), similar to (19), describes a vibrating membrane, but with different elastic stiffness and a non-uniform distribution of mass across the transformed region. In contrast, for the model of a flexural plate, equation (9), after the transformation (5), does not preserve the physical interpretation, i.e. it is no longer the equation of free vibrations of a plate. This suggests that the problem in hand is very different from the model of a cloak for a membrane. It presents an additional challenge to identify the physical configuration consistent with the new equation (9). This issue is discussed in the next section.

The cloaking transformation does not produce an orthotropic inhomogeneous plate
We make a direct comparison between the transformed equation (9), and the equation for an inhomogeneous orthotropic plate. Firstly, we note that the moments θ M M , r and θ M r satisfy the partial differential equation: For an orthotropic inhomogeneous plate, where rigidities and Poissonʼs coefficients vary radially, equation (22) has the form ⎛ Direct comparison of (23) and (9) shows that the fourth-order terms agree in the equation of the orthotropic plate and the transformed equation within the cloaking region. However, a discrepancy occurs in other lower-order terms, and hence the transformed equation (9) does not represent a classical orthotropic Kirchhoff plate. Additional physical constraints are needed to complete the model. This will be achieved through the approximation discussed in the next section.

The cloaking approximation
Equation (23) can be rewritten after the substitution of flexural rigidity coefficients as in (8): Direct comparison with equation (7) shows the discrepancy in the third-order derivative terms in addition to that in the lower-order terms. The difference between the left-hand sides of equations (7) and (24) is  It is apparent that all coefficients in the above equation have the form f r R r ( ) j 1 , with f j being smooth functions when > r R 1 , and these coefficients are small when R r 1 is considered as a small parameter, in particular, when the penetration depth for the incident wave into the cloaking region is small.
In the approximation implemented here, we chose the parameters of the cloak in such a way that the interior diameter of the cloaking ring is sufficiently small compared to the diameter of the whole cloaking region and compared to the wavelength of the incident wave, i.e. the following non-dimensional quantities are small: where β is defined by (13). The material outside the cloak remains unaffected by the transformation, whereas the interior material represents a radially orthotropic plate in the framework of the approximation described here (see equations (24), (25)). Numerical simulations below show the efficiency of our concept, which is also in agreement with the experimental evidence published in [18].

Numerical illustration
The notion of an approximate cloak, introduced above, is used here in the numerical illustrations. This approximation is valid for a certain choice of geometrical parameters and frequency values. Numerical simulations are produced for an elastic, isotropic Kirchhoff plate which contains a radially orthotropic, inhomogeneous cloaking layer (the flexural displacement of the plate on the interior of cloaking layer is governed by (24)). Without loss of generality, the incident field is represented by a flexural plane wave propagating horizontally. Perfectly matched layers (PML) are used on the exterior boundary of the computational domain. PML conditions are 'absorbing' boundary conditions simulating a non-reflective exterior contour. The parameters used for the numerical simulations are shown in table 1. The exterior of the cloak corresponds to a homogeneous, isotropic plate, whereas the interior of the cloak is an inhomogeneous, radially orthotropic plate. The numerical simulations were produced using Comsol Multiphysics ® (see the appendix for more details on the numerical implementation).
In figure 1, we consider the case of interior and exterior radii for the cloaking region chosen as R 1 = 0.   It is also expected that the approximation is frequency sensitive, and the properties of the approximate cloak may also change with the variation of the thickness of the cloak. This is illustrated in figure 2. In part (a) of that figure, the simulation corresponds to the case of a higher frequency (ω = 200), and the cloaked obstacle shows a non-suppressed shadow. Similarly, in part (b) of figure 2 we have non-suppressed shadow for a different reason. Although the frequency of the incident wave remains the same as in figure 1, the size of the obstacle has increased and the interior radius of the cloak is twice as large as the case in figure 1(b). Consequently, in both diagrams shown in figure 2 the cloaking has been affected. 6. Alternative approach: plate subjected to in-plane forces and pre-stress In this section we show that, by choosing a different normalization, it is possible to give a physical interpretation of the transformed plate equations as a Kirchhoff plate subjected to inplane forces in addition to the usual flexural behavior. This can also lead to a broadband perfect cloak. Here, we extend the cartesian formulation given recently in Colquitt et al [19] to the cylindrical cloak configuration. In particular, equations (7) and (10) are normalized in the following way: Cloaking has been affected in these two cases.  Then, introducing the following definition for the rigidity and inertial parameters: These quantities are constrained to satisfy the in-plane balance equation The final form for the transformed equation is leading to a consistent physical interpretation. In [19], for a different cloak geometry, we have shown that such a pre-stressed elastic system leads to broadband cloaking. We note the resemblance of the above computations in figures 1 and 2 with those in [18], which shows results from an experimental study of a structured cloak and flexural waves. In [18] the cloaking approximation is shown to be frequency sensitive, so that cloaking does not occur for frequencies above a certain threshold. From (1)-(3), it is clear that four independent elastic parameters D r , θ D , θ D r , and ν r are required to characterize a radially orthotropic plate (also see the classic papers [30,31]). Only Youngʼs moduli E r and θ E appear to be given in [18].
The inertial properties are defined by the mass density which is also required when implementing the cloak, both in computations and experiments. Different normalization of the mass density can be applied; in particular, the mass density used in [18,32] was constant. Here we have defined all of the required parameters and explained how they fit into the configuration approximating the flexural cloak. We have also given the range of validity of such an approximation.

Asymptotic derivation of the transformed plate equation from the equations of elasticity
In this section, the transformed equation of motion for the Kirchhoff plate (10) is deduced directly from the transformed equations of motion of three-dimensional linear elasticity. An asymptotic model is implemented in order to obtain the lower-dimensional plate model from the analysis of a thin three-dimensional solid. It was shown in earlier works [4,7,16] that the transformed equations of elasticity are subject to the choice of gauge. In particular, the resulting material may lack the minor symmetries in the constitutive equations. This does not occur in the case of flexural plates, as demonstrated below.

Transformed equations of elasticity
The Navier equations  In (32), χ Ω = × −h h [ 2,2], with  Ω ⊆ 2 , λ and μ are the Lamé moduli, ρ is the mass density of the medium and zero body forces are assumed. Field equations (32) are accompanied by homogeneous Neumann boundary conditions on the upper and lower external surfaces = ± Z h 2: . Now, we introduce a geometric transformation Accordingly, the Jacobi matrix F (in cylindrical coordinates) and the Jacobian J are given by Then, the Navier equations (32) transform into ⎡

Asymptotic model
In order to obtain the Kirchhoff plate model directly from the transformed equations of elasticity (36) and (37) an asymptotic procedure for elliptic operators in thin domains is developed [33,34]. We introduce the scaled spatial variable ξ ϵ = z , ϵ ≪ 1, and we also assume that the transverse displacement component depends on the scaled time variable ϵ = T t. Then, we consider the following asymptotic approximation for the displacement vector u after substitution of the equations of motion (36) and boundary conditions (37). The solvability condition for ξ W (0) constitutes a well-posed problem for the transverse displacement field ξ v (0) describing the flexural behavior of a thin plate. The solvability conditions for W r (0) and θ W (0) constitute a well-posed problem for the in-plane displacement field  r (0) and  θ (0) describing the behavior of a thin shell. Here interest is in the description of the plate model and we restrict attention to the asymptotic procedure for ξ v (0) which will be indicated by v for ease of notation. After the introduction of the scaled variables ξ and T, the equations of motion (36) and boundary conditions (37) can be expressed in the form .
In equation (39)    and the differential operator  R 2 1 is defined in equation (11). For equation (40) where v does not depend on ξ and the solvability conditions are automatically satisfied.
To the next order, the field equations and boundary conditions Note that v (1) , and v (2) and v (3) in the following, are complemented by the normalization condition of zero average along the thickness. Next, the field equations and boundary conditions accompanied by the boundary conditions (2) and the corresponding solution is Finally, the vector function W (0) satisfies the equation 2 (0) 2 together with the boundary conditions In particular, ξ W (0) solves the problem ⎛ where E is Youngʼs modulus, ν Poissonʼs ratio and H and T have been replaced by ϵ h and ϵ t, respectively.

Conclusions
There are different ways of reducing the shadow generated by a scatterer. In particular, an elementary example where a 'heavy' inclusion is surrounded by a 'lighter' isotropic coating was discussed in [35]. That model requires the average mass density of the inclusion and coating together to be the same as the mass density of the ambient matrix. Such examples have been known for more than a century (see, for example, [36]). It is important to mention that a combination of a heavy inclusion and a lighter coating cannot be associated with an 'invisibility cloak', but instead can be used to reduce the monopole source term in the asymptotics at infinity.
In [18] the use of micro-structured material for cloaking represented a substantial advance. That work has demonstrated that cloaking of a flexural wave is possible, although such a cloaking approximation is frequency dependent. In the present paper, we have provided a full theoretical background for such an approximation and have also discussed the range of its applicability.
Furthermore, by referring to pre-stressed elastic plates, we have resolved a long-standing problem of creating an exact cloak for flexural waves. For the cloaking region obtained as a result of a 'push-out' radially symmetric transformation, we have identified a full set of parameters, including pre-stress and in-plane body forces. In the case when pre-stress and body forces are not included in the model, an approximation of the cloak has been developed for ≪ R R 1 1 2 and within the frequency range when β ≪ R 1 1 . The illustrative numerical computations show excellent agreement with the prediction of the theoretical model and the existing experimental results.
The transformed equations of three-dimensional vector elasticity were analysed asymptotically for a thin solid. The resulting lower-dimensional model agrees fully with the outcome of the direct application of the radial 'push-out' transformation to the equation of motion of a Kirchhoff plate. It is also noted that 'transformed' material in three-dimensional elasticity has non-symmetric constitutive relations, as outlined in [16], but the lowerdimensional model for the plate does not have such a feature. The physical nature of the reduced model is fully explained, with the introduction of pre-stress and in-plane body forces, which have been identified in explicit closed form. Implementation of the proposed model could lead to a new generation of lightweight and highly-efficient structured shields and filtering devices.
Kichhoff-Love plate equation may be simulated in the finite element code using the Mindlin-Reissner Comsol package. For a detailed comparison of the dynamics of Kirchhoff-Love and Mindlin plates we refer to [39].
A benchmark example is considered here to illustrate the choice of parameters in the Comsol computational model. This involves a radially symmetric plate loaded by a point force that is applied at the centre of the plate. This problem has a closed-form solution, and as expected, the Comsol numerical scheme shows excellent agreement in the framework of the Kirchhoff model.
For the purpose of the numerical simulations presented in section 5.1, the shear modulus was chosen as ≈ G 10 10 and = × − h 1 10 3 , while all other parameters were chosen as unity. In order to verify this approach, Comsolʼs Mindlin plate model was used to compute a static verification model. In particular, Greenʼs function for a homogeneous radially orthotropic circular plate of radius R 2 with clamped boundaries was considered. This problem was considered for Kirchhoff-Love plates in [40], section 82 (see also [31]), and subject to correction of typographical mistakes, has the following analytical solution