Numerical Assessment of the Performance of Elastic Cloaks for Transient Flexural Waves

A relevant application of transformation elastodynamics has shown that flexural waves in a Kirchhoff-Love plate can be diverted and channeled to cloak a region of the ambient space. To achieve the goal, an orthotropic meta-structural plate should be employed. However, the corresponding mathematical transformation leads to the presence of an unwanted strong compressive prestress, likely beyond the buckling threshold of the structure, with a set of in-plane body forces to warrant equilibrium. In addition, the plate must possess, at the same time, high bending stiffnesses, but a null twisting rigidity. With the aim of estimating the performance of cloaks modelled with approximate parameters, an in-house finite element code, based on a subparametric technique, is implemented to deal with the cloaking of transient waves in orthotropic thin plates. The tool allows us to explore the sensitivity of specific stiffness parameters that may be difficult to match in a real cloak design. In addition, the finite element code is extended to investigate a meta-plate interacting with a Winkler foundation, to confirm how the subgrade modulus should transform in the cloak region.


State of the Art and Research Challenges
The control of elastic waves to cloak a region of the ambient space has been shown achievable by transformation elastodynamics, which provides the mechanical properties of the material surrounding the region (Milton et al., 2006;Brun et al., 2009;Norris and Shuvalov, 2011). A relevant application of this broad area is the control of transverse waves in plates, for which solutions have been proposed mainly based on two approaches: a "passive" one, where the features of the cloak are achieved by a given microstructure (Farhat et al., 2009a;Farhat et al., 2009b;Brun et al., 2014;Colquitt et al., 2014;Climente et al., 2016;Zareei and Alam, 2017;Liu and Zhu, 2019;Misseroni et al., 2019) and an "active" one, in which tunable quantities depending on the actual mechanical input are employed for the same goal (see e.g., Futhazar et al., 2015;O'Neill et al., 2015;Chen et al., 2016;Ning et al., 2020). Experiments have been also proposed to validate some of the previous theoretical proposals (Stenger et al., 2012;Misseroni et al., 2016;Darabi et al., 2018;Misseroni et al., 2019). Most of these investigations show that this technique can be feasible and likely to be employed in centimetre-size real-life systems with wave frequencies in the order of a few kHz.
Within the set of attempts based on the "passive" approach, Colquitt et al. (2014) proposed a transformation to conceive a Kirchhoff-Love plate theory for cloaking flexural waves. The resulting governing equations for the thin plate in the transformed domain involve the presence of variable (in space) bending and torsional stiffnesses, and density, in addition to the presence of in-plane body forces and prestresses in the plate. In the same paper, the general framework was specialised to the case of a square cloak composed of four trapezoidal elements (see Figure 1) embedded in an isotropic, homogeneous domain. For this geometry, a set of relationships was established to provide explicit expressions of the quantities concerned in each part of the cloak. The broadband effectiveness of the meta-material plate cloak was then assessed by means of numerical tests.
The features of the suggested meta-structure can be summarised as follows: (1) the cloak is locally an orthotropic thin plate with principal directions varying point-to-point. These directions obey the geometric symmetries of the domain (see, e.g., trapezoid displayed in Figure 1B, where the principal directions are sketched with thin lines). In the same figure, the increase of the out-of-alignment of the local principal system with the axis of symmetry towards the diagonals of the cloak can be noticed; (2) by investigating how the stiffnesses of the plate in the orthotropic principal system depend on the position and thinking to construct the plate with a homogeneous material, it turns out that the thickness of the plate must assure an increase in bending stiffness along the direction parallel to the inner boundary of the cloak moving away the axis of symmetry of each trapezoid and a decrease of the bending stiffness along the "radial" direction moving towards the centre of the cloak. These two requirements are apparently in contradiction, but they must be addressed in an effective design; (3) the twisting stiffness in the principal system of orthotropy vanishes at each point of the domain; (4) the mass density of the cloak is not constant and varies with the Jacobian of the transformation; in particular, this quantity decreases by approaching the inner boundary of the cloak; (5) in-plane body forces and prestresses in the interior of the domain as well as forces per unit length applied along the diagonals of the cloak are necessary to warrant equilibrium. The predicted prestress is compressive in a large part of the domain with values that are likely to exceed the buckling threshold for the thin structure.
The set of listed properties demonstrates that the design and engineering of an effective cloak for flexural waves based on transformation elastodynamics of Kirchhoff-Love plates is an exceptional challenge. Actually, the presence of severe compressive prestresses leads to the conclusion that it is impossible to construct a stable, ideal cloak. Therefore, an approach based on carefully considered approximations are required to achieve the goal. In the experimental test conducted by Misseroni et al. (2016), for instance, the design of the approximated cloak has focused mainly on the match of bending stiffnesses along the directions of orthotropy.
Moving from the points raised in the just edited list, the goal of the paper is four-fold: (1) the first aim is to design and implement a fully open in-house FE code, based on a subparametric technique, to simulate transient wave propagation in locally orthotropic Kirchhoff-Love plates; (2) the numerical tool is then exploited to assess the performance of meta-structural square cloaks designed by assuming reasonable approximations (prestress free, local twisting stiffness as limited as possible, cross sections of the plate ensuring the required local bending stiffnesses, non-constant mass density); (3) the same tool is employed to investigate the role of certain stiffness parameters, i.e., twisting and coupling bending stiffnesses; (4) in addition, a non-secondary purpose of the paper is to propose an extension of the theory to encompass interaction of the cloak with an elastic substrate modelled with a Winkler foundation and assess its effectiveness.

Background
The definition of the flexural cloak is based on a transformation, say χ, that maps a two-dimensional subdomain C 0 ⊂ K occupied by a homogeneous, isotropic Kirchhoff-Love plate, to the domain C where the meta-plate is to be constructed, i.e., χ k : The map is such that the transformed equation in the latter is still the governing equation for flexural waves supported by a Kirchhoff-Love plate. In our case, C 0 is a square deprived of a tiny hole H 0 in the centre ( Figure 1A) that is mapped to a collection of four trapezoids C (i) (i 1, . . . , 4) ( Figure 1B), such that C ∪ C (i) . The length of one side of the external boundary of both C 0 and C is equal to 2L whereas that of H 0 is equal to 2εa, where ε is a small parameter. 2 In both configurations, the exterior domain K∖C 0 K∖C is occupied by the initial, homogeneous structure that possesses constant thickness H, mass density per unit volume ρ 0 and bending stiffness D 0 EH 3 /[12(1 − ] 2 )], where E is the Young's modulus of the material and ν its Poisson's ratio. When free-standing, the governing equation of the plate under harmonic vibrations takes the form where w(x 0 K , t) is the transverse displacement, ∇ 4 0 is the biharmonic operator in C 0 and σ 0 ρ 0 H is the mass density per unit area. For the relevant case of a plate resting on a homogeneous substrate modelled as a bed of independent linear springs (Winkler model), a term k 0 W is to be added on the left-hand side of Eq. 1, where k 0 W is the subgrade coefficient. We emphasise that in Eq. 1, D 0 and σ 0 are independent of the position.
The function χ k (x 0 K ) is assumed to be invertible, a requisite that is met if its evaluation in the inner boundary of the cloak is non-singular as in the case under study, After having recalled that uppercase indices refer to points belonging to C 0 and that a comma indicates partial differentiation, the gradient of the transformation x i,K , its Jacobian J(x p ) det x i,K and the tensor 3 g ij x p x i,P x j,P J are instrumental in defining the properties of the cloak as shown by Colquitt et al. (2014).
In order to have a complete overview on the established theory of meta-plate cloaks, we remind that the local governing equation of a Kirchhoff-Love plate takes the form (Timoshenko and Woinowsky-Krieger, 1959;Lekhnitskii et al., 1968) where σ is the mass density per unit area, m ij represents the set of moments per unit length whilst n ij and s i denote membrane forces and in-plane body forces per unit length, respectively, in equilibrium through fulfillment of the condition n ij,j + s i 0.
The next step is to substitute the constitutive equations that connect m ij to the generalised curvature w ,kl as m ij −D ijkl (x p )w ,kl , where D ijkl is the stiffness tensor which may depend on the position. The substitution in Eq. 2 yields an equation used by Colquitt et al. (2014) to identify the different terms corresponding to the transformation and confirm that the flexural displacement of a Kirchhoff-Love plate under an arbitrary coordinate mapping may be interpreted as a generalised plate. The identification of terms shows that the stiffness tensor of the transformed plate corresponds to the in-plane body forces to 1 Uppercase and lowercase indices, both ranging between 1 and 2, are referred to reference and transformed domains, respectively, whereas a 0 in either super-or sub-script position indicates that the quantity concerned is evaluated in the reference configuration. 2 ϵ plays the role of regularisation parameter for the construction of a near cloak (Kohn et al., 2008) in which the material properties at the inner boundary of the cloak are not singular. 3 The summation over the repeated index is implied throughout the paper.
Frontiers in Materials | www.frontiersin.org November 2020 | Volume 7 | Article 603667 s l x p D 0 g rs Jg ql,q ,r ,s , and the transformed density per unit area turns out to be σ(x p ) ρ 0 H/J. When the Winkler foundation is involved, the subgrade coefficient transforms similarly to the density, i.e., a quantity that is now dependent on the position. Focusing now on the square cloak transformation sketched in Figure 1 and, in particular, on trapezoid C (1) , and {e 1 , e 2 } is the orthonormal basis associated with both coordinate systems Ox 0 1 x 0 2 and Ox 1 x 2 ; a is the half-length of the side of the square hole to be cloaked and b L − a is the "cloak thickness". This information defines the six independent components of the stiffness tensor D ijkl , the prestress n ij and the in-plane body forces s i in the four subdomains. For trapezoid C (1) , all these quantities are reported in Appendix A. For the other parts of the cloak, they can be easily inferred by taking into account the relevant symmetries. Note that α is dimensionless For the same trapezoid C (1) , to reach the orthotropic principal directions at a point P, the angle θ(x 1 , x 2 ) (anticlockwise in the domain x 2 > 0, see Figure 1B) is now introduced. In particular, this angle identifies the principal axis x (respectively y) moving from x 1 (respectively x 2 ). Its value is where In the new local system Pxy, the constitutive relationships assume the form where the new stiffness moduli contained therein are provided in terms of D ijkl in Appendix A.
Interestingly, as shown by Colquitt et al. (2014), the twisting stiffness D xyxy brought about by the transformation vanishes in all points of the transformed domain. This is a particular feature of the cloak under investigation that must be properly addressed in a meaningful design. It is worth recalling that in an isotropic plate with compact cross section, D xyxy D 0 (1 − ]).
The cloak is subjected to traction-free boundary conditions in the inner free boundary and to continuity conditions at the interface between itself and the outer homogeneous plate.

NUMERICAL MODEL AND IMPLEMENTATION IN A IN-HOUSE FINITE ELEMENT CODE
With the aim of performing numerical simulations via the Finite Element Method (FEM) of the transient dynamic behaviour of a square cloak, a weak form of the governing equations is now obtained. The formulation can be derived by multiplying Eq. 4 by the test function v(x 1 , x 2 ) and using the Galerkin method (Zienkiewicz and Taylor, 2005;Reddy, 1993). After the application of the divergence theorem to the generic subdomain Ω of the entire plate and use of the equilibrium Eq. 3, the final expression becomes where the term involving k W is now taken into account. In expression Eq. 11, (·) ,n and (·) ,ŝ are the directional derivatives along the normal unit vector n and tangent unit vector s of the boundary Γ of Ω; m n , m nŝ are the bending and twisting moments per unit length along Γ, respectively, whereas n n , n nŝ describe the in-plane internal actions along the normal and tangent unit vectors; t n is the specific shear force on the boundary and H 2 0 (Ω) is the suitable Sobolev space. It should be noted that in the standard procedure leading to the weak form Eq. 11, the terms D ijkl,i and D ijkl,ij naturally disappear by integrating by part. The analysed plate is composed of various subdomains with different stiffness properties represented by the terms D ijkl : the uncloaked region K∖C possesses a homogeneous stiffness D 0 whereas the cloak is described by the parameters D ijkl given in Eq. 5.
Geometric and analytical approximations are performed to achieve the discretised formulation by dividing the entire domain into a finite subset of elements, i.e., Ω ≈ ∪ n e e 1 Ω e , where n e is the number of the elements, and substituting unknown and test functions with the linear combinations (12) respectively, where N k (x 1 , x 2 , t), k {p, q}, are shape functions defined on the spatial domain Ω e and indices p, q range within a Frontiers in Materials | www.frontiersin.org November 2020 | Volume 7 | Article 603667 suitable interval as specified later. For each Finite Element, the final discretised formulation can be written as where n p and n q are the number of the shape functions that discretise w and v, respectively, and all the terms related to the inplane prestress have been neglected. A parametric representation is implemented in the code: the shape functions are defined in the local coordinate system Oξη of the master element Ω e {(ξ, η)4R 2 : −1 ≤ ξ ≤ 1,−1 ≤ η ≤ 1}, which is mapped into each element of the mesh using appropriate coordinate transformations. Since an isoparametric representation is not possible for thin-plate elements unless a previous distortion of the mesh (see for example Petera and Pittman, 1994), a subparametric representation has been adopted, so the bilinear shape functions that are used for the transformation of the master element are different from the hermitian ones used for the analytical approximation of the unknown.
The structured mesh employed in the numerical analyses is composed of quadrangular elements with four nodes in the corner of each element whose generalised displacements are ,ξη , k 1, 2, 3, 4, leading to (n p n q ) 16 degrees of freedom for the element (Imbert, 1979;Petyt, 1990;Reddy, 1993;Zienkiewicz and Taylor, 2005). Finally, after substitution and integration of the aforementioned shape functions in Eq. 1 and the assembly operation over all the Finite Elements of the mesh, the set of algebraic equations used for the spatial approximation of this problem can be written as where M and K are the mass matrix and the stiffness matrix, respectively; u(t) is the vector of the degrees of freedom whereas f (t) is the vector of generalised external loads. The approximation in the time domain is performed with the single-step time-integration algorithm Generalized-α method (Chung and Hulbert, 1993;Chung and Hulbert, 1994) which is an evolution of the well-known Newmark method that is frequently used for transient dynamical analysis in mechanics (see Imbert, 1979;Hughes, 1987;Petyt, 1990;Zienkiewicz and Taylor, 2005). According to this method, Eq. 15 is written as where a generic damping matrix C comes into play and the timedomain discretisation T {0, t 1 , t 2 , . . . , t n , . . . , t N T max }, t n+1 t n + Δt n , is performed with respect to the time-domain T [0, T max ]. The Taylor's expansions of the solution u(t n+1 ) ≡ u n+1 and its first time derivative _ u(t n+1 ) ≡ _ u n+1 are given by where θ 1 and θ 2 are positive real parameters that define the integration scheme. Eq. 18 are substituted into Eq. 16 where α f and α m are two additional parameters that will be specified later. The following implicit scheme can be used to achieve the solution € u n+1 , namely where Δt 2 n 2 . Equations 18 and 20 are employed to calculate the complete solution at each time step t n .
The Generalized-α method is completely described by the four parameters α m , α f , θ 1 , θ 2 that define the stability and the convergence properties of the time-integration scheme. They can be written as a function of the asymptotic spectral radius ρ ∞ ∈ [0, 1] which is related to the numerical damping in the highfrequency limit as It is worth mentioning that Eq. 15 can be profitably employed to perform time-harmonic simulations once a harmonic form for u(t) and f (t) is introduced.
The mass matrix can generally be computed in different ways: the consistent mass matrix M C is obtained directly from the weak Eq. 13, so the inertia contribution is assigned to each node via the shape function; the lumped mass matrix M L is a diagonal matrix obtained condensing the mass of each element in its nodes via equilibrium. In the ensuing simulations we considered a mixed mass matrix M (1 − α l )M C + α l M L , where α l ∈ [0, 1] [see Hughes (1987) and Zienkiewicz and Taylor (2005)].
Frontiers in Materials | www.frontiersin.org November 2020 | Volume 7 | Article 603667 5 We implemented this numerical model in a Finite Element code, written using the software Matlab. The program deals with anisotropic and inhomogeneous thin elastic plates with different possible shapes of the domain and various kind of boundary and loading conditions. Using this code, it is also possible to perform static and modal analyses.

NUMERICAL SIMULATIONS AND PERFORMANCE OF THE CLOAK
Using the code described above, we studied the transient propagation of flexural waves in a simply-supported square plate in which a square hole is cloaked with respect to out-ofplane vibration generated by a point force varying sinusoidally. With reference to Figure 2, the homogeneous plate in C 0 is a steel one (E 210 GPa, ] 0.3, ρ 0 7, 800 kg/m) with thickness H 1 mm and bending stiffness D 0 19.23 Nm. The outer side of the cloak is such that L 0.75 m while a 0.1923 m, b 0.5577 m and εa 0.025 m. 4 When a Winkler-like soil comes into play in the analysis, the stiffness of the elastic substratum is represented by k 0 W 0.5 · 10 6 N/m, qualitatively corresponding to that of a silicone material (Bigoni et al., 2008;Gei, 2008).
The boundary conditions on the outer perimeter are those of a simply supported plate, namely null transversal displacement w| zK 0, null rotation w ,ŝ | zK 0 and null bending moment m n | zK 0. The distance between the centre of the plate and the centre of the cloaked square hole is equal to 2.5 m. A harmonic transversal point force F(t) F cos(ωt) is applied at the centre of the plate where we assumed F 100 N and ω varying from 50 to 600 rad/s. The information of the point load is contained in the vector load f (t) in the right-hand-side of Eqs. 15 and 16.
The mesh that we used for the transient analysis is composed of (n e ) 83,340 quadrangular elements with 335,864 total degrees of freedom; the elements are rectangular and regular outside the area of the cloak whereas the mesh is refined and composed of trapezoidal elements inside the area of the cloak. Each trapezoid of the cloak C (i) is composed of 1,080 Finite Elements.
We performed a transient dynamical analysis with the Generalized-α method (see above) for which we chose ρ ∞ 1 to neglect numerical damping, therefore the method reduced to a second order and unconditionally stable algorithm with α m α f θ 1 θ 2 1/2, which is equivalent to the constant-average acceleration method (Reddy, 1993). This choice of the parameters of the numerical algorithm leads to a stable scheme for any time interval, so we considered a time domain T [0, 200 ms] with Δt n 1 ms. We assumed homogeneous initial condition, so u(t 0) u 0 0 and _ u(t 0) _ u 0 0. The initial value of € u 0 (necessary in the first step) has been calculated assuming the initial condition on u 0 and _ u 0 and solving Eq. 16 for t 0. In the time-domain integration we neglected the contribution of the damping matrix, i.e., C 0, and we assumed α l 0.5. The performance of the cloak is now assessed with reference to the following case studies: (1) a cloak possessing all properties of the Theoretical Transformation Without Prestress (TTWP) is analysed first. As recalled in the Introduction, the prestress is ruled out as it is not realistic to assume the distribution of initial stress n ij and in-plane body forces foreseen by the transformation. For this cloak, however, the local twisting stiffness D xyxy vanishes in the whole domain C; (2) a cloak possessing all properties of the TTWP-but with a different amount of torsional stiffness D xyxy to evaluate the influence of this parameter-is investigated. This case is significant as every constructed meta-plate cloak would be equipped with a certain amount of twisting stiffness; (3) same as 2), but with the goal of estimating the role of the coupling bending stiffness D xxyy ; (4) a cloak with the properties predicted by the TTWP is then studied, but resting on a substrate modelled as a Winkler foundation whose position-dependent subgrade modulus obeys relationship Eq. 7.
For all case studies, the displacement maps and the other relevant quantities that are illustrated are computed at a time t compatible with a propagation of the wavefront located downline of the region of the cloak, but not significantly disturbed by secondary waves reflected at the boundary of the plate.

Performance of the Cloak
The displacement map of a TTWP cloak computed at t 155 ms for ω 300 rad/s is displayed in Figure 3B, to be compared with that reported in Figure 3A for a plate with an uncloaked hole. In the latter diagram, a wake is evident just on the right-hand side of the hole as well as a perturbed wave pattern on the left-hand side due to reflection caused by the hole. Those two features have been eliminated by the cloak that is able to regularise the wavefront emerging from the device into the expected circular pattern and minimise back scattering. A cloak following the transformation with prestress (Eqs. 5 and 6) would have produced a circular wavefront matching that of a homogeneous plate as presented by Colquitt et al. (2014).
In order to appraise more in detail the quality of the cloaking, Figure 4 reports the out-of-plane displacement w for the four sections sketched in Figure 2; the first three sections are transverse with respect to the propagation of the wave, at a distance from the centre of the hole of 1.25 m (Section no. 1), 1.75 m (Section no. 2) and 3 m (Section no. 3); the fourth one runs longitudinally from a point at 0.25 m just outside the righthand boundary of the cloak to the external side of the plate. In all plots, the red line, which describes the displacements of the TTWP cloak, follows quite well the dashed line representing the response of the square, hole-less, homogeneous plate K. Figure 4D shows that the approximate cloak is able to reproduce quite satisfactorily the phase of the wave of the homogeneous plate/perfect cloak at points in the region just after the cloaked region, to confirm the remark on Figure 3B that the TTWP cloak is able to reconstruct the circular wavefront.
In order to measure quantitatively the quality of the cloak we adopt the index ) that is computed along the relevant section ς that is described by the co-ordinate r. Perfection corresponds to Q 0. In Eq. 22, w FE are the displacements computed by the numerical model for the non-homogeneous case under investigation and w Hom corresponds to the solution for the homogeneous plate. With reference to the four sections displayed in Figure 4, it turns out that the quality indices are those reported in Table 1, showing that the cloaked solutions are in average 5.95 times more effective than the uncloaked ones, therefore a good result is in any case obtained. The effect of an amount of twisting stiffness on the displacement map of a cloak TTWP is displayed in Figure 5 where several increasing values of the parameter D xyxy are analysed. The plot on the far left is the one reported in Figure 3B. The detrimental effect of this parameter is readily observed. A quality factor can also be computed for the displayed cases. For Section no. 2, Q increases from 0.1241 (second panel) to 1.348 (far right panel) showing a constant worsening in the response of the meta-plate. In Appendix B, a strategy to minimise the twisting stiffness of a fibre-composite micro-structured plate is illustrated. A similar analysis can be reiterated for the coupling stiffness D xxyy . This term plays an unexpected important role as revealed by the simulations reported in Figure 6 where TTWP metaplates are implemented, but possessing a D xxyy whose point-wise value corresponds to 90, 50, and 10% of the theoretical one. Unexpected because, on the one hand, in a rational microstructured plate design where the principal bending stiffnesses D xxxx and D yyyy are matched first, the coupling term simply results as an outcome of the choices made previously; on the other hand, the mathematical transformation generally predicts values of this parameter that are remarkably larger than those that can be reasonably reachable for an orthotropic plate (see Appendix B, where the list of theoretical stiffnesses are reported for two selected cases). The three panels displayed in the figure clearly demonstrate the relevance of the coupling bending stiffness in assuring an accurate behaviour of the cloak.  The last picture of this subsection (Figure 7) refers to a TTWP cloak subjected to point-loads pulsating at different ω, to show that the approximate model is behaving well for a quite broad range of frequencies. The time of computation t differs from one picture to the other in order to compare a similar displacement pattern. For all displayed cases, the meta-plate is able to reconstruct satisfactorily the circular wavefront.

Flexural Cloaks on Winkler Substrates
In real applications, cloaks are likely grounded to a substrate which also provides the support of the object to hide. The cloaking features of a plate resting on a Winkler foundation is analysed in Figure 8 for ω 300 rad/s at t 145 ms. In part b), the TTWP cloak rests on a bed of springs whose subgrade modulus is constant (k W k 0 W ). This assumption clearly worsens the performance of the device with respect to the uncloaked hole case ( Figure 8A). Conversely, the position dependent stiffness k W (x p ), whose value in C follows Eq. 7, guarantees a satisfactory behaviour of the device, as it can be noted that the wavefront reconstructs correctly after the transformed domain. Similarly to the simply-supported plate, solutions along two sections (i.e. no. 2 and no. 4 in Figure 2A) are FIGURE 5 | Displacement map w (in m) computed numerically at t 155 ms for ω 300 rad/s carried out for increasing twisting stiffness D xyxy . The detrimental effect on the cloaking performance of an increase in this parameter is clearly evident. 0.3304, respectively. Again, the behaviour of the perfect cloak, both in phase and amplitude, is captured remarkably well by the approximate one with position dependent subgrade modulus. This shows that a correct design of a flexural cloak resting on a substrate requires the presence of a variable subgrade modulus that follows the provision of the theoretical transformation.
As a final remark of this section, one must note that the employed geometrical, elastic and loading conditions correspond to a chosen prototype geometry, but the scale of the dimensions, and the load intensity and frequency can be changed to obtain an

CONCLUDING REMARKS
The engineering of a meta-structural cloak for elastic flexural waves based on transformation elastodynamics is an exceptional challenge which theoretically would require the implementation of unfeasible compressive prestresses and in-plane body forces to warrant equilibrium. Therefore, reasonable assumptions should be adopted that are however grounded on the following main points: 1) the meta-plate has a locally orthotropic response with 2) an almost vanishing twisting stiffness, 3) spatially-varying bending stiffnesses and 4) density. In addition, the structure invariably would interact with a substrate.
With the aim of dealing with a fully open simulation tool, a FE code is developed and implemented with the specific purpose of studying transient wave propagation in locally orthotropic Kirchhoff-Love plates. A subparametric technique is adopted for spatial discretization, whereas the approximation in the time domain is performed with the single-step time-integration algorithm Generalized-α method.
The performance of a prototype meta-plate square cloak based on the approximate assumption of null prestress is then assessed parametrically by focusing specially on the role of the twisting stiffness and the coupling bending stiffness. The reason is that while the local principal bending rigidities can be matched quite easily in a microstructured plate, for those two parameters the process is much more difficult.
The simulations show that for an effective cloaking response, the twisting stiffness should be lower than 10% of the bending one for the homogeneous plate, while for the coupling bending stiffness a departure of 10% from the theoretical value already worsens with some evidence the performance of the meta-plate.
A second contribution of this work consists in the extension of the general theory of thin-plate cloaks to comprise interaction of the meta-structure with an elastic substrate via Winklerfoundation model. The conclusion is that the subgrade modulus transforms similarly to the mass density of the plate. The numerical simulations confirm this finding and clearly show that a constant modulus beneath the cloak jeopardizes dramatically the functionality of the structure.
The FE tool can be further expanded to embed inelastic responses of plate and substrate, thus enabling simulations of approximate supported meta-plates subjected to transient waves arising from seismic shocks.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation

AUTHOR CONTRIBUTIONS
DV and MG conceived the idea. MR and DV developed and implemented the numerical tool. MR performed the simulations. MR, DV, and MG wrote the paper.  1) the function c c B (x 1 , x 2 ) is computed through Eq. 25 in which D xxxx (x 1 , x 2 ) replaces the current entry on the left-hand side of the equation. The resulting expression is a cubic in the unknown. For the typical involved parameter D xxxx (x 1 , x 2 ), there is always a solution in the range 0 < c c B ≤ 0.8; 2) Equation 25 can now be employed to compute H y or alternatively c t B (x 1 , x 2 )) in which D yyyy (x 1 , x 2 ) appears on the left-hand side of the equation. In a practical case, the usual choice is to set c t B 0.8 and then H y is to be calculated; 3) the teeth density n/H should be selected so that the twisting stiffness can be estimated via Eq. 25.
Eventually, the coupling term D clk xxyy can be computed through Eq 25.
The outcome of the design is shown in Tables 2 and 3 with reference to two transformations whose parameters are (both with L 0.75 m): 1) a 0.1 m, b 0.65 m and εa 0.025 m, and 2) a 0.1923 m, b 0.5577 m and εa 0.025 m (i.e., the one studied in the numerical simulations). It can be noted that in all control points listed in the tables, the value of D clk xxyy is orders of magnitude smaller than the theoretical one.
FIGURE 10 | Local micro-structure of the meta-plate for cloaking flexural waves. The yellow (resp. pink) cross section provides the bending stiffness D xxxx (resp. D yyyy ). Axes x and y correspond to the local principal directions of orthotropy. The teeth density is n/H.   0.1923, 0.1923) 10.71 -11.37 --- The calculations are for n 5. a At this point, the principal system of orthotropy is aligned with Ox1x2. b At this point, the principal directions of orthotropy are rotated of an angle θ C 0.282 rad (16 + 171) with respect to Ox1x2.