Skip to main content

ORIGINAL RESEARCH article

Front. Mater., 27 November 2020
Sec. Mechanics of Materials
This article is part of the Research Topic Cutting-Edge Metastructures: Micro-Architected and Active Metamaterials View all 8 articles

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

Marco RossiMarco Rossi1Daniele VeberDaniele Veber1Massimiliano Gei,
&#x;Massimiliano Gei1,2*
  • 1Department of Civil, Environmental and Mechanical Engineering, University of Trento, Trento, Italy
  • 2School of Engineering, Cardiff University, Cardiff, United Kingdom

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.

Introduction

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.

FIGURE 1
www.frontiersin.org

FIGURE 1. Sketch of a square meta-structural cloak studied by Colquitt et al. (2014): (A) initial, undeformed domain C0 where the plate is isotropic and homogeneous; (B) transformed domain C where the four trapezoids composing the cloak are numbered. The principal directions of orthotropy are represented with thin lines (the local axes at a generic point P are labelled x and y as depicted). The perimeter of trapezoid C(1) is marked with a dashed red line.

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 C0K occupied by a homogeneous, isotropic Kirchhoff-Love plate, to the domain C where the meta-plate is to be constructed, i.e., χk:C0C,xK0xk=χk(xK0).1 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, C0 is a square deprived of a tiny hole H0 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 C0 and C is equal to 2L whereas that of H0 is equal to 2εa, where ε is a small parameter.2

In both configurations, the exterior domain K∖C0=KC is occupied by the initial, homogeneous structure that possesses constant thickness H, mass density per unit volume ρ0 and bending stiffness D0=EH3/[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

D004w+σ0w¨=0,(1)

where w(xK0,t) is the transverse displacement, 04 is the biharmonic operator in C0 and σ0=ρ0H 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 kW0 is to be added on the left-hand side of Eq. 1, where kW0 is the subgrade coefficient. We emphasise that in Eq. 1, D0 and σ0 are independent of the position.

The function χk(xK0) 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 C0 and that a comma indicates partial differentiation, the gradient of the transformation xi,K, its Jacobian J(xp)=detxi,K and the tensor3

gij(xp)=xi,Pxj,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)

mij,ij+nijw,ijsiw,i=σw¨,(2)

where σ is the mass density per unit area, mij represents the set of moments per unit length whilst nij and si denote membrane forces and in-plane body forces per unit length, respectively, in equilibrium through fulfillment of the condition

nij,j+si=0.(3)

The next step is to substitute the constitutive equations that connect mij to the generalised curvature w,kl as mij=Dijkl(xp)w,kl, where Dijkl is the stiffness tensor which may depend on the position. The substitution in Eq. 2 yields

Dijklw,ijkl+2Dijkl,iw,jkl+(Dijkl,ijnkl)w,kl+slw,l=σw¨,(4)

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

Dijkl(xp)=D0Jgijgkl,(5)

the in-plane body forces to

sl(xp)=D0[grs(Jgql,q),r],s,(6)

and the transformed density per unit area turns out to be σ(xp)=ρ0H/J.

When the Winkler foundation is involved, the subgrade coefficient transforms similarly to the density, i.e.,

kW(xp)=kW0/J,(7)

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), χ(1)(xK0eK)=(αx10+c)e1+(αx20+cx20/x10)e2,

xi,K(1)=(α0x2αcx1(cx1)x1αx1c),

and J(1)=x1α2/(x1c), where

α=b/L1εa/L, c=aεa1εa/L,

and {e1,e2} is the orthonormal basis associated with both co-ordinate systems Ox10x20 and Ox1x2; a is the half-length of the side of the square hole to be cloaked and b=La is the “cloak thickness”. This information defines the six independent components of the stiffness tensor Dijkl, the prestress nij and the in-plane body forces si 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 while [c]=[L].

For the same trapezoid C(1), to reach the orthotropic principal directions at a point P, the angle θ(x1,x2) (anticlockwise in the domain x2>0, see Figure 1B) is now introduced. In particular, this angle identifies the principal axis x (respectively y) moving from x1 (respectively x2). Its value is

θ=14arccos[QQ+2(D1111D2222)(D2221D1112)],(8)

where

Q=(D1111+D22222D11224D1212)(D1112+D2221).(9)

In the new local system Pxy, the constitutive relationships assume the form

mxx=Dxxxxw,xxDxxyyw,yy, myy=Dxxyyw,xxDyyyyw,yy, mxy=Dxyxyw,xy,(10)

where the new stiffness moduli contained therein are provided in terms of Dijkl in Appendix A.

Interestingly, as shown by Colquitt et al. (2014), the twisting stiffness Dxyxy 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, Dxyxy=D0(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(x1,x2) and using the Galerkin method (Zienkiewicz and Taylor, 2005; Reddy, 1993). After the application of the divergence theorem to the generic sub-domain Ω of the entire plate and use of the equilibrium Eq. 3, the final expression becomes

Ωv,ijDijklw,kldΩΩv,lnklw,kdΩ+Ωσvw¨dΩ+ΩkWvwdΩΓ{mn^v,n^v(tn^+mn^s^,s^)}ds+Γ{nn^s^vw,n^+nn^vw,s^}ds=0,  vH02(Ω),(11)

where the term involving kW is now taken into account. In expression Eq. 11, (),n^ and (),s^ are the directional derivatives along the normal unit vector n^ and tangent unit vector s^ of the boundary Γ of Ω; mn^, mn^s^ are the bending and twisting moments per unit length along Γ, respectively, whereas nn^, nn^s^ describe the in-plane internal actions along the normal and tangent unit vectors; tn^ is the specific shear force on the boundary and H02(Ω) is the suitable Sobolev space. It should be noted that in the standard procedure leading to the weak form Eq. 11, the terms Dijkl,i and Dijkl,ij naturally disappear by integrating by part. The analysed plate is composed of various subdomains with different stiffness properties represented by the terms Dijkl: the uncloaked region KC possesses a homogeneous stiffness D0 whereas the cloak is described by the parameters Dijkl 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., Ωe=1neΩe, where ne is the number of the elements, and substituting unknown and test functions with the linear combinations

w(x1,x2,t)=Np(x1,x2,t)up(t),v(x1,x2,t)=Nq(x1,x2,t)vq,(12)

respectively, where Nk(x1,x2,t), k={p,q}, are shape functions defined on the spatial domain Ωe and indices p,q range within a suitable interval as specified later. For each Finite Element, the final discretised formulation can be written as

p=1npq=1nqvq{ΩeNq,ijDijklNp,klupdΩ+ΩeσNqNpu¨pdΩ+ΩekWNqNpupdΩΓe[mn^Nq,n^Nq(tn^+mn^s^,s^)]ds}=0,vq,(13)

where np and nq are the number of the shape functions that discretise w and v, respectively, and all the terms related to the in-plane 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={(ξ,η)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

{w(k)w,ξ(k)w,η(k)w,ξη(k)},k=1,2,3,4,(14)

leading to (np=nq=) 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

Mu¨(t)+Ku(t)=f(t),(15)

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

Mu¨(tn)+Cu˙(tn)+Ku(tn)=f(tn),(16)

where a generic damping matrix C comes into play and the time-domain discretisation

T={0,t1,t2,,tn,,tN=Tmax},tn+1=tn+Δtn,(17)

is performed with respect to the time-domain T=[0,Tmax].

The Taylor’s expansions of the solution u(tn+1)un+1 and its first time derivative u˙(tn+1)u˙n+1 are given by

un+1=un+Δtnu˙n+(1θ2)Δtn22u¨n+θ2Δtn22u¨n+1,u˙n+1=u˙n+(1θ1)Δtnu¨n+θ1Δtnu¨n+1,(18)

where θ1 and θ2 are positive real parameters that define the integration scheme. Eq. 18 are substituted into Eq. 16 calculated in the intermediate time value tn+1αf=(1αf)tn+1+αftn, thus obtaining

M[(1αm)u¨n+1+αmu¨n]+C[(1αf)u˙n+1+αfu˙n]+K[(1αf)un+1+αfun]=fn+1αf,(19)

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

u¨n+1=A1[fn+1αfKun(C+K(1αf)Δtn)u˙n(Mαm+C(1αf)(1θ1)Δtn+K(1αf)(1θ2)Δtn22)u¨n],(20)

where A=[M(1αm)+C(1αf)θ1Δtn+K(1αf)θ2Δtn22]. Equations 18 and 20 are employed to calculate the complete solution at each time step tn.

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 high-frequency limit as

αm=2ρ1ρ+1, αf=ρρ+1, θ1=12αm+αf, θ2=12(1αm+αf)2.(21)

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 MC 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 ML 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)MC+αlML, where αl[0,1] [see Hughes (1987) and Zienkiewicz and Taylor (2005)].

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-of-plane vibration generated by a point force varying sinusoidally. With reference to Figure 2, the homogeneous plate in C0 is a steel one (E=210 GPa, ν=0.3, ρ0=7,800 kg/m) with thickness H=1 mm and bending stiffness D0=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 kW0=0.5106 N/m, qualitatively corresponding to that of a silicone material (Bigoni et al., 2008; Gei, 2008).

FIGURE 2
www.frontiersin.org

FIGURE 2. A) Geometry of the studied plate where it is visible the point where the time-varying force is enforced (dimensions in m) and the sections considered in Figure 4; (B) detail of the mesh, the external boundary of the cloak is marked with a dashed red line.

The boundary conditions on the outer perimeter are those of a simply supported plate, namely null transversal displacement w|K=0, null rotation w,s^|K=0 and null bending moment mn^|K=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 (ne=) 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 Δtn=1 ms. We assumed homogeneous initial condition, so u(t=0)=u0=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 u0 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 nij and in-plane body forces foreseen by the transformation. For this cloak, however, the local twisting stiffness Dxyxy vanishes in the whole domain C;

(2) a cloak possessing all properties of the TTWP–but with a different amount of torsional stiffness Dxyxy 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 Dxxyy;

(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).

FIGURE 3
www.frontiersin.org

FIGURE 3. Displacement map w (in m) computed numerically at t=155 ms for ω=300 rad/s: (A) uncloaked square hole; (B) performance of a cloak obtained through the Theoretical Transformation Without Prestress (TTWP). The square-shaped white contour sketched in (B) indicates the external boundary of the cloak.

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 right-hand 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.

FIGURE 4
www.frontiersin.org

FIGURE 4. Displacement map w (in m) computed at t=155 ms for ω=300 rad/s along the four Sections sketched in Figure 2. Comparison between solutions for a cloaked (TTWP–red line) and an uncloaked (blue line) hole. The dashed line represents the response of a homogeneous plate without the hole.

In order to measure quantitatively the quality of the cloak we adopt the index (Colquitt et al., 2014)

Q=ς|wFE(r)wHom(r)|2drς|wHom(r)|2dr,(22)

that is computed along the relevant section ς that is described by the co-ordinate r. Perfection corresponds to Q=0. In Eq. 22, wFE are the displacements computed by the numerical model for the non-homogeneous case under investigation and wHom 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.

TABLE 1
www.frontiersin.org

TABLE 1. Quality index for the four sections reported in Figure 4.

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 Dxyxy 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.

FIGURE 5
www.frontiersin.org

FIGURE 5. Displacement map w (in m) computed numerically at t=155 ms for ω=300 rad/s carried out for increasing twisting stiffness Dxyxy. The detrimental effect on the cloaking performance of an increase in this parameter is clearly evident.

A similar analysis can be reiterated for the coupling stiffness Dxxyy. This term plays an unexpected important role as revealed by the simulations reported in Figure 6 where TTWP meta-plates are implemented, but possessing a Dxxyy 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 Dxxxx and Dyyyy 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.

FIGURE 6
www.frontiersin.org

FIGURE 6. Displacement map w (in m) computed numerically at t=140 ms for ω=300 rad/s carried out for decreasing coupling bending stiffness Dxxyy. DxxyyTTWP indicates the value of the variable predicted by the theoretical transformation. The damaging effect on the cloaking performance of a decrease in this parameter is clearly evident.

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.

FIGURE 7
www.frontiersin.org

FIGURE 7. Displacement maps w (in m) computed numerically at frequencies 100, 200, 300 and 400 rad/s (at different times t in order to compare a similar displacement pattern) for a TTWP cloak.

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 (kW=kW0). This assumption clearly worsens the performance of the device with respect to the uncloaked hole case (Figure 8A). Conversely, the position dependent stiffness kW(xp), 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 studied in Figure 9. The quality index for the two sections are Q2cl=0.0134, Q2uncl=0.1719 and Q4cl=0.0164, Q4uncl=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.

FIGURE 8
www.frontiersin.org

FIGURE 8. Displacement maps w (in m) computed numerically at t=145 ms for ω=300 rad/s. (A) Uncloaked square hole where the plate rests on a Winkler foundation with kW0=0.5106 N/m; (B) hole cloaked by a plate obtained through the TTWP resting on a Winkler foundation with a constant kW=kW0; (C) same as (B), but with a position-dependent subgrade module whose value follows Eq. 7. The square-shaped white contour indicates the external boundary of the cloak.

FIGURE 9
www.frontiersin.org

FIGURE 9. Displacement map w (in m) computed at t=145 ms for ω=300 rad/s along Sections no. 2 (A) and no. 4 (B) in Figure 2 for a plate resting on a Winkler foundation with kW0=0.5106 N/m. Comparison between solutions for a cloaked hole (TTWP with position-dependent subgrade module (see Eq. 7 and Figure 8) –red line) and an uncloaked one (constant subgrade module–blue line). The dashed line represents the response of a homogeneous plate without the hole.

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 identical behaviour on a different scaled structure subjected to different load conditions.

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 Winkler-foundation 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.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

MR acknowledges Italian Ministry of Education, University and Research (MIUR) for supporting their research under grant PRIN no. 2015LYYXA8. MG acknowledges support from EU-H2020-Marie Curie-Sklodowska Action-COFUND grant SIRCIW (agreement no. 663830).

Appendix A: Material Parameters and Pre-Stress for the Cloak

The six independent flexural stiffnesses for trapezoid C(1) of the cloak are

D1111(1)=α2(1cx1)D0,D2222(1)=α2(c2x22+x14)2(x1c)3x15D0,D2211(1)=α2(c2x22+x14)(x1c)x13D0,
D1212(1)=α2c2x22(x1c)x13D0,D1112(1)=α2cx2x12D0,D2212(1)=α2cx2(c2x22+x14)(x1c)2x14D0.

The remaining components can be deduced from the major and minor symmetries of Dijkl. The membrane forces and in-plane body forces are

n11(1)=2α2cx12(x1c)D0,  n22(1)=2α2c(x14+8cx22x13c2x22)x14(x1c)3D0,
n12(1)=2α2cx2(3x12c)(x1c)2x13D0,s1(1)=0,s2(1)=24α2cx2(x1c)3x12D0.

Stiffnesses transform locally as follows:

Dxxxx=D1111cos4θ+2(D1122+2D1212)sin2θcos2θ+D2222sin4θ+2(D1112cos2θ+D2221sin2θ)sin2θ,
Dyyyy=D1111sin4θ+2(D1122+2D1212)sin2θcos2θ+D2222cos4θ2(D1112sin2θ+D2221cos2θ)sin2θ,
D1=D1122+(D1111+D22222(D1122+2D1212))sin2θcos2θ+(D2221D1112)cos2θsin2θ,
Dxyxy=D1212+(D1111+D22222(D1122+2D1212))sin2θ cos2θ+(D2221D1112)cos2θsin2θ.

Appendix B: An Example of a Micro-Structured Plate Cloak

A design based on a micro-structured meta-plate is here proposed with the aim of matching locally the principal bending stiffnesses Dxxxx and Dyyyy, and limiting the twisting stiffness Dxyxy. The exercise is conducted by assuming a high-performance fibre-reinforced material whose epoxy matrix (shortened as “ep”; material parameters: Eep=3.4 GPa, νep=0.3, ρep=1,200 kg/m) is stiffened by long boron fibres (shortened as “B”; material parameters: EB=380 GPa, νB=0.13, ρB=2,600 kg/m) (Kaw, 2006; Mroz, 2011). The meta-plate is such that the cross section orthogonal to axis x is compact (and corresponds to the “core”, indicated with “c”) with its thickness matching the height H and fibres aligned along axis x. That orthogonal to axis y features a set of tiny rectangular appendices, in the number of n over a length equal to H, on both the outer sides called “teeth” (shortened as “t”). The total height of core and teeth is Hy (Figure 10). The reason of selecting rectangular appendices lies in the stringent requirement of limiting the twisting stiffness. The slenderness of the teeth may lead to instability on the side of the plate in compression, however this issue is not further addressed here. Fibres are here aligned along axis y only in the teeth. Along both directions, fibres are locally measured out (their volume fractions are spatially-varying design variables that are however capped at 0.8) to allow the relevant cross section to reach the needed bending stiffness.

In the spirit of the Kirchhoff-Love theory for plates, the core material is subjected to a plane-stress state, therefore the linear elastic constitutive equations read

σxx=Ec11νxycνyxcεxx+νxycEc21νxycνyxcεyy,  σyy=Ec21νxycνyxcεyy+νxycEc21νxycνyxcεxx,
τxy=Gcγxy,(23)

whereas the teeth undergo a uniaxial stress (i.e., σyyt=Etεyy). All constitutive parameters of the composite in both core and teeth follow the rule of mixtures applied to composite materials,5 however more sophisticated models [e.g., Halpin-Tsai model, Halpin and Kardos (1976)] can be adopted for their estimation.

The effective stiffnesses Dijklclk appearing in the cloak moment/curvature constitutive equations (cf. Eq. 10)

mxx=Dxxxxclkw,xxDxxyyclkw,yy,  myy=Dyyyyclkw,yyDxxyyclkw,xx,(24)
mxy=Dxyxyclkw,xy

can be calculated by applying the standard methodology of integrating locally stresses across the plate thickness. The integration leads to

Dxxxxclk=D0REc1Ec2,  Dxxyyclk=D0Rνxyc,  Dxyxyclk=GcH36+GtHyH6(Hn)2,
Dyyyyclk=D0[R+(1ν2)EtE[(HyH)31]],(25)

where R=(Ec2/E)(1ν2)/(1νxycνyxc).

With reference to the square transformation displayed in Figure 1B, the properties of the two cross sections of the micro-structured plate can be defined as follows:

1) the function cBc(x1,x2) is computed through Eq. 25 in which Dxxxx(x1,x2) 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 Dxxxx(x1,x2), there is always a solution in the range 0<cBc0.8;

2) Equation 25 can now be employed to compute Hy or alternatively cBt(x1,x2)) in which Dyyyy(x1,x2) appears on the left-hand side of the equation. In a practical case, the usual choice is to set cBt=0.8 and then Hy 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 Dxxyyclk 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 Dxxyyclk is orders of magnitude smaller than the theoretical one.

FIGURE 10
www.frontiersin.org

FIGURE 10. Local micro-structure of the meta-plate for cloaking flexural waves. The yellow (resp. pink) cross section provides the bending stiffness Dxxxx (resp. Dyyyy). Axes x and y correspond to the local principal directions of orthotropy. The teeth density is n/H.

TABLE 2
www.frontiersin.org

TABLE 2. Microstructural parameters of the boron/epoxy cloak for the transformation whose parameters are a=0.1 m, b=0.65 m, εa=0.025 m, L=0.75 m.

TABLE 3
www.frontiersin.org

TABLE 3. Microstructural parameters of the boron/epoxy cloak for the transformation whose parameters are a=0.1923 m, b=0.5577 m, εa=0.025 m, L=0.75 m.

Footnotes

1Uppercase 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.

3The summation over the repeated index is implied throughout the paper.

4With this set of parameters, it turns out that for the trapezoid C(1) the stiffnesses for the four representative points sketched in Figure 1B are displayed in Table 3. The minimum value D1111 is achieved along the whole inner boundary whereas the maximum D2222 is at point D.

5In particular, Ec1, νxyc follow the weighted average ()ξ=()BcBξ+()epcepξ, whereas Ec2, Gc, and Gt obey the harmonic average 1/()ξ=cBξ/()B+cepξ/()ep, with ξ{c,t}; νyxc=νxycEc2/Ec1.

References

Bigoni, D., Gei, M., and Movchan, A. B. (2008). Dynamics of a prestressed stiff layer on an elastic half space: filtering and band gap characteristics of periodic structural models derived from long-wave asymptotics. J. Mech. Phys. Solid. 56, 2494–2520. doi:10.1016/j.jmps.2008.02.007

CrossRef Full Text | Google Scholar

Brun, M., Colquitt, D. J., Jones, I. S., Movchan, A. B., and Movchan, N. V. (2014). Transformation cloaking and radial approximations for flexural waves in elastic plates. New J. Phys. 16, 093020. doi:10.1088/1367-2630/16/9/093020

CrossRef Full Text | Google Scholar

Brun, M., Guenneau, S., and Movchan, A. B. (2009). Achieving control of in-plane elastic waves. Appl. Phys. Lett. 94, 061903. doi:10.1063/1.3068491

CrossRef Full Text | Google Scholar

Chen, Y., Hu, J., and Huang, G. (2016). A design of active elastic metamaterials for control of flexural waves using the transformation method. J. Intell. Mater. Syst. Struct. 27, 1337–1347. doi:10.1177/1045389x15590273

CrossRef Full Text | Google Scholar

Chung, J., and Hulbert, G. M. (1993). A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-α method. J. Appl. Mech. 60, 371–375. doi:10.1115/1.2900803

CrossRef Full Text | Google Scholar

Chung, J., and Hulbert, G. M. (1994). A family of single-step houbolt time integration algorithms for structural dynamics. Comput. Methods Appl. Mech. Eng. 118, 1–11. doi:10.1016/0045-7825(94)90103-1

CrossRef Full Text | Google Scholar

Climente, A., Torrent, D., and Sánchez-Dehesa, J. (2016). Analysis of flexural wave cloaks. AIP Adv. 6, 121704. doi:10.1063/1.4968611

CrossRef Full Text | Google Scholar

Colquitt, D. J., Brun, M., Gei, M., Movchan, A. B., Movchan, N. V., and Jones, I. S. (2014). Transformation elastodynamics and cloaking for flexural waves. J. Mech. Phys. Solid. 72, 131–143. doi:10.1016/j.jmps.2014.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Darabi, A., Zareei, A., Alam, M.-R., and Leamy, M. J. (2018). Experimental demonstration of an ultrabroadband nonlinear cloak for flexural waves. Phys. Rev. Lett. 121, 174301. doi:10.1103/physrevlett.121.174301

PubMed Abstract | CrossRef Full Text | Google Scholar

Farhat, M., Guenneau, S., and Enoch, S. (2009a). Ultrabroadband elastic cloaking in thin plates. Phys. Rev. Lett. 103, 024301. doi:10.1103/PhysRevLett.103.024301

PubMed Abstract |CrossRef Full Text | Google Scholar

Farhat, M., Guenneau, S., Enoch, S., and Movchan, A. B. (2009b). Cloaking bending waves propagating in thin elastic plates. Phys. Rev. B 79, 033102. doi:10.1103/PhysRevB.79.033102

CrossRef Full Text | Google Scholar

Futhazar, G., Parnell, W. J., and Norris, A. N. (2015). Active cloaking of flexural waves in thin plates. J. Sound Vib. 356, 1–19. doi:10.1016/j.jsv.2015.06.023

CrossRef Full Text | Google Scholar

Gei, M. (2008). Elastic waves guided by a material interface. Eur. J. Mech. Solid. 27, 328–345. doi:10.1016/j.euromechsol.2007.10.002

CrossRef Full Text | Google Scholar

Halpin, J. C., and Kardos, J. L. (1976). Halpin-tsai equations: a review. Polym. Eng. Sci. 16, 344–352. doi:10.1002/pen.760160512

CrossRef Full Text | Google Scholar

Hughes, T. (1987). The finite element method - linear static and dynamic finite element analysis. Upper Saddle River, NJ: Prentice-Hall.

Google Scholar

Imbert, F. (1979). Analyse des structures par elements finis. Paris, France: Cepadues editions.

Google Scholar

Kaw, A. (2006). Mechanics of composite materials. Boca Raton, FL: Taylor & Francis Group.

Google Scholar

Kohn, R. V., Shen, H., Vogelius, M. S., and Weinstein, M. I. (2008). Cloaking via change of variables in electric impedance tomography. Inverse Probl. 24, 015016. doi:10.1088/0266-5611/24/1/015016

CrossRef Full Text | Google Scholar

Lekhnitskii, S., Tsai, S., and Cheron, T. (1968). Anisotropic plates. New York, NY: Gordon and Breach.

Google Scholar

Liu, M., and Zhu, W. D. (2019). Nonlinear transformation-based broadband cloaking for flexural waves in elastic thin plates. J. Sound Vib. 445, 270–287. doi:10.1016/j.jsv.2018.12.025

CrossRef Full Text | Google Scholar

Milton, G. W., Briane, M., and Willis, J. R. (2006). On cloaking for elasticity and physical equations with a transformation invariant form. New J. Phys. 8, 248. doi:10.1088/1367-2630/8/10/248

PubMed Abstract | CrossRef Full Text | Google Scholar

Misseroni, D., Colquitt, D. J., Movchan, A. B., Movchan, N. V., and Jones, I. S. (2016). Cymatics for the cloaking of flexural vibrations in a structures plate. Sci. Rep. 6, 23929. doi:10.1038/srep23929

PubMed Abstract | CrossRef Full Text | Google Scholar

Misseroni, D., Movchan, A. B., and Bigoni, D. (2019). Omnidirectional flexural invisibility of multiple interacting voids in vibrating elastic plates. Proc. R. Soc. A 475, 2229. doi:10.1098/rspa.2019.0283

PubMed Abstract | CrossRef Full Text | Google Scholar

Mroz, A. (2011). Stability analysis of a plane, rectangular, boron-epoxy laminated plate basing on strength properties determined by different methods. Mech. Mech. Eng. 15, 161–181.

Google Scholar

Ning, L., Wang, Y.-Z., and Wang, Y.-S. (2020). Active control cloak of the elastic wave metamaterial. Int. J. Solid Struct. 202, 126–135. doi:10.1016/j.ijsolstr.2020.06.009

CrossRef Full Text | Google Scholar

Norris, A. N., and Shuvalov, A. L. (2011). Elastic cloaking theory. Wave Motion. 48, 525–538. doi:10.1016/j.wavemoti.2011.03.002

CrossRef Full Text | Google Scholar

O’Neill, J., Selsil, O., McPhedran, R. C., Movchan, A. B., and Movchan, N. V. (2015). Active cloaking of inclusions for flexural waves in thin elastic plates. Q. J. Mech. Appl. Math. 68, 263–288. doi:10.1093/qjmam/hbv007

CrossRef Full Text | Google Scholar

Petera, J., and Pittman, J. F. T. (1994). Isoparametric hermite elements. Int. J. Numer. Meth. Engng. 37, 3489–3519. doi:10.1002/nme.1620372006

CrossRef Full Text | Google Scholar

Petyt, M. (1990). An introduction to the finite element method. Cambridge, UK: Cambridge University Press.

Google Scholar

Reddy, J. (1993). An introduction to the finite element method. 2nd Edn. New York, NY: McGraw-Hill.

Google Scholar

Stenger, N., Wilhelm, M., and Wegener, M. (2012). Experiments on elastic cloaking in thin plates. Phys. Rev. Lett. 108, 014301. doi:10.1103/PhysRevLett.108.014301

PubMed Abstract | CrossRef Full Text | Google Scholar

Timoshenko, S., and Woinowsky-Krieger, S. (1959). Theory of plates and shells. New York, NY: McGraw-Hill.

Google Scholar

Zareei, A., and Alam, M.-R. (2017). Broadband cloaking of flexural waves. Phys. Rev. E 95, 063002. doi:10.1103/PhysRevE.95.063002

PubMed Abstract | CrossRef Full Text | Google Scholar

Zienkiewicz, O., and Taylor, R. (2005). The finite element method. 6th Edn. Amsterdam, Netherlands: Elsevier.

Google Scholar

Keywords: meta-materials, transformation elastodynamics, invisibility cloak, wave propagation, finite element method

Citation: Rossi M, Veber D and Gei M (2020) Numerical Assessment of the Performance of Elastic Cloaks for Transient Flexural Waves. Front. Mater. 7:603667. doi: 10.3389/fmats.2020.603667

Received: 07 September 2020; Accepted: 12 October 2020;
Published: 27 November 2020.

Edited by:

Andrea Bacigalupo, University of Genoa, Italy

Reviewed by:

Francesca Fantoni, University of Brescia, Italy
Marco Lepidi, University of Genoa, Italy

Copyright © 2020 Rossi, Veber and Gei. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Massimiliano Gei, massimiliano.gei@dia.units.it

†Current address: Department of Engineering and Architecture, University of Trieste, Trieste, Italy.

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.