The paper presents an approach to modelling a novel elastic metamaterial structure that possesses non-trivial dispersion features. A system of spinners has been embedded into a two-dimensional periodic lattice system. The analysis of the motion of the spinners is used to derive an expression for a ‘chiral term’ in the equations describing the dynamics of the lattice. Dispersion of elastic waves is shown to possess innovative filtering and polarization properties induced by the vortex-type nature of the structured media. The related effective behaviour in a continuous medium is implemented to build a shielding cloak around an obstacle. Analytical work is accompanied by numerical illustrations.
Interaction of electromagnetic waves with photonic crystal structures has been the subject of rapid development in the last decade. In particular, cloaking of electromagnetic waves was analysed both experimentally and theoretically; one of the earliest investigations in this direction was reported by Schurig et al. (2006). Further influential study in this area was published by Li & Pendry (2008). The concept of cloaking was also addressed for other types of governing equations, such as surface water waves (Farhat et al. 2008) and vibration of plates (Farhat et al. 2009). It is noted, however, that the equations of vector elasticity pose a considerable challenge compared with Maxwell's system or the bi-harmonic operator for plate elastic structures: in simple terms, this is linked to the tensor structure of governing equations, which involve the notion of shear. As a consequence, there are two types of waves propagating with different speeds coupled to each other via the boundary conditions. A recent commentary on this subject was published by McPhedran & Movchan (2012). The modelling of elastic metamaterial systems, which possess cloaking properties, remains a challenge and this paper presents a constructive algorithm addressing this issue.
Dynamic response of vector elastic lattice systems, interacting with waves, has been the subject of classical investigations (e.g. Slepyan 1981; Kunin 1982), in addition to the more recent publications of Slepyan (2002), Slepyan & Ayzenberg-Stepanenko (2002) and Colquitt et al. (2011). We note substantial differences between scalar problems involving vibration of systems of harmonic springs and vector problems of elasticity, referring to elastic rods and beams, which connect a system of finite solids or point masses. The notion of shear stress and shear strain is absent in scalar lattice systems and there is no such analogue in models of electromagnetism. The misconception of predictability of dynamic properties of vector elastic lattices is sometimes based on an intuitive extrapolation of results available for the scalar systems or for problems of electromagnetism in which there is only one wave speed. A thorough analysis is the only correct way forward and, as illustrated in Colquitt et al. (2011), some results for micro-polar elastic lattice structures do not follow from the simpler physical models or intuitive assumptions.
Phononic band-gap structures have been analysed in recent years for periodic arrays of voids or finite inclusions in the continuum matrix by Movchan et al. (1997), Poulton et al. (2000), Zalipaev et al. (2002), Platts et al. (2002) and Lin & Huang (2011). Platts et al. (2003a,b) studied transmission problems for arrays of elastic structured stacks including a comparative analysis of the filtering properties of elastic waves in doubly periodic media and the transmission properties for the corresponding singly periodic stack structure. In particular, a review of the Rayleigh multiple methods was given by McPhedran et al. (2001).
Furthermore, a highly cited paper (Milton et al. 2006) addresses a specially designed elastic coating, with anisotropic inertia properties, which may act as an ‘invisibility cloak’, i.e. it routes a wave around an obstacle. Brun et al. (2009) show an alternative in the design of an invisibility cloak which incorporates a micro-polar composite.
In this paper, we address propagation of elastic waves, their dispersion properties, and interaction with defects in a two-dimensional structured medium supplemented with micro-rotations created by spinning masses embedded into the lattice system. The analysis of gyroscopic motion of an individual mass is incorporated into the system of conservation of linear and angular momenta within the lattice systems. Furthermore, elastic Bloch–Floquet waves are considered in such a system. It is shown that the presence of spinning masses ‘stiffens’ the overall system, and it also leads to a special design of chiral media, which are characterized by ‘handedness’ in their micro-structures. These media may possess stop-band properties in addition to shielding (re-routing) with respect to elastic waves. In turn, we consider several types of continuous chiral systems to design coatings around finite defects in such a way that a dynamic signature of the coated defect may be enhanced or suppressed compared with the case when the coatings are absent.
2. A vortex-type lattice system
A general algorithm was established for analysing the spectral properties of a periodic lattice in Martinsson & Movchan (2003), and lattice systems with built-in dynamic rotational interactions were analysed in Colquitt et al. (2011). The new model proposed in this paper incorporates the vortex type of action at the junction regions of the lattice system. We also derive equations containing the vortex terms and explain their physical nature.
(a) Governing equations
We consider a triangular periodic lattice with the elementary cell shown in figure 1, which is defined in terms of the basis vectors, 2.1In this case, we have two junctions per elementary cell of the periodic structure, and it is assumed that the point masses placed at these junctions are m1 and m2. The stiffnesses of elastic truss-like links, each of length l, are assumed to be equal to c, and the mass density of the elastic links is assumed to be negligibly small.
We also use the three unit vectors a(j),j=1,2,3, characterizing the directions of trusses within the elementary cell: The multi-index n=(n1,n2) with integer entries is used to characterize the position of the elementary cell within the periodic structure, so that the position vector of the mass mj (j=1,2) is defined as 2.2where the vectors t(1) and t(2) are given in (2.1) and the displacement u(n,j) eiωt is time-harmonic, of radian frequency ω. For convenience, we also use the multi-indices e1=(1,0) and e2=(0,1).
The design of the lattice system considered here includes gyroscopic spinners attached to every junction. The axis of each spinner is perpendicular to the plane of the two-dimensional lattice, and its angular velocity is related (see (2.10)) to the radian frequency ω of the oscillation of the main lattice system, as outlined in §2b. A small in-plane displacement of a junction leads to a change in orientation of the spinner axis and hence generates a moment creating a ‘vortex-type’ effect.
The equations of motion for the masses m1 and m2 within the elementary cell n can be formed as in Martinsson & Movchan (2003). In our case, they include gyroscopic terms where the gyroscopic force is orthogonal to the displacement 2.3and 2.4where R is the rotation matrix and κ1, κ2 are the spinner constants.
(b) Evaluation of the spinner constants
It is assumed that each junction within the lattice is connected to a spinner, whose axis is perpendicular to the plane of the lattice, as shown in figure 1. By considering one of the junctions, together with the spinner attached, we evaluate the spinner constants in the form κj=αj (i,j=1,2,), where αj are real coefficients that have the dimension of mass.
Let ψ, ϕ and θ be the angles of spin, precession and nutation with respect to the vertical axis Oz, respectively (e.g. Goldstein et al. 2000, p. 210). An individual spinner and the angles ψ, ϕ and θ are shown in figure 1c. Then, assuming that gravity is absent, the equations of motion of the spinner, for the case of the constant spin rate (), are 2.5 2.6 and 2.7where Mx,My,Mz are the moments about the x,y and z axes, and I0=Ixx=Iyy,I=Izz are the moments of inertia.
We consider a small amplitude time-harmonic motion (owing to the motion of the lattice attached to the spinner), resulting in the nutation angle and constant spin rate .
Assuming that there are no imposed moments about the x and the y axes (i.e. Mx=My=0) and the precession rate is constant (), we deduce the connection between the spin rate Ω and the radian frequency ω of the imposed nutation. The precession rate is given as 2.8According to the equation of motion (2.5), we have 2.9The compatibility condition for equations (2.8) and (2.9) has the form 2.10Direct substitution into (2.7) gives the induced moment about the z- axis Let h be the characteristic length of the spinner, as shown in figure 1c. Then the magnitude of the in-plane displacement in the lattice junction is Taking into account that the moment about the z-axis is Mz=Fhθ, the ‘rotational force’ F is given by which indicates the presence of a phase shift in the gyroscopic force compared with the time-harmonic displacement. Following this derivation, the spinner constants α1 and α2 in equations (2.3) and (2.4) are
(c) Elastic Bloch waves
Consider the case when m1=m2 and α1=α2, i.e. the monatomic lattice. Then equations (2.3) and (2.4) are equivalent, subject to a translation of an elementary cell of the periodic system, which includes only a single junction. This configuration is a good example for observing the influence of the vortex-type interaction within the lattice on the dispersion properties of the elastic waves. The amplitudes U1,U2 of the time-harmonic displacement, of radian frequency ω, within the periodic vortex-type lattice satisfy the system of equations (2.3), together with the Bloch–Floquet conditions set within the elementary cell of the periodic system, as follows: 2.11where k is the Bloch vector (k1,k2), and matrix T has the form Formula (2.11) refers to the case of a monatomic lattice where the elementary cell is smaller than the case of a bi-atomic system (compare with (2.2)). These waves are dispersive, and the corresponding dispersion equation (see Martinsson & Movchan 2003) is 2.12where the mass matrix M=diag[m,m] and Σ represents the chiral term associated with the presence of spinners attached to the junction points of the lattice system. Assuming that each spinner is characterized by the ‘spinner constant’ α, we represent the matrix Σ as 2.13
Additionally, in equation (2.12), C(k) is the stiffness matrix for the monatomic lattice of masses m, with spring connectors of stiffness c. It is given by where 2.14
We note that the dispersion equation is bi-quadratic with respect to ω, and it has the form 2.15
Since c>0, both the trace trC and the determinant detC are positive for all k1 and k2 in the elementary cell of the reciprocal lattice, except at the origin, where both of them are zero. Hence, there are two important regimes. For m2>α2, there are two dispersion surfaces. As the factor (m2−α2) decreases towards a critical point where m2=α2, the upper dispersion surface increases without limit and the lower surface decreases. At this critical point, the upper surface is infinite and the lower surface has the degenerate value of (m−1 detC/trC)1/2. For m2<α2, the lower dispersion surface only remains. Representative results will now be given for each of the regimes m2>α2 (subcritical) and m2<α2 (supercritical).
(d) Dispersion properties
Firstly, figure 2 gives dispersion surfaces and the corresponding two-dimensional diagrams with the slowness contours for the uniform lattice without spinners (α=0). Two conical surfaces are clearly visible in the low-frequency range. These surfaces are shown to evolve as the spinners are brought into the system, now with chiral terms in the equations of motion.
(i) Monatomic lattice of the vortex type
The diagrams for the vortex type of monatomic lattice with a small spinner constant1 α=0.5 and mass m=1 are given in figure 3. It is apparent that the lower acoustic surface, dominated by shear waves, does not change much as a result of the vortex interaction. On the contrary, the upper acoustic surface, dominated by pressure waves, extends further into a higher range of frequencies and has a pronounced conical shape for a vortex type of lattice of small α. The latter implies that the homogenization range for Bloch–Floquet waves corresponding to this dispersion surface is also extended into a wider range of frequencies, as shown in figure 3.
For the other regime, the dispersion diagram illustrating such a situation for the case of α=2 and m=1 is given in figure 4. Only a single dispersion surface is seen for the low-frequency range.
The influence of the vortex type of interaction appears to be significant for the dispersion properties of Bloch–Floquet waves within the entire admissible frequency range. The illustration is given in figure 5, where we assume k1=k2=k and plot ω as a function of the spinner constant α and the Bloch parameter k. In this case, the lower dispersion surface is represented by the diagram in figure 5a and is defined for all real values of α, whereas the upper dispersion surface disappears in the supercritical regime when |α|≥m. The diagram characterizing the upper dispersion surface and showing ω as a function of α and k is presented in figure 5b for 0≤α<m, where m=1. It is noted that, for a fixed value of k, the frequency ω on the upper dispersion surface increases with an increase in the spinner constant α (figure 5b), whereas on the lower dispersion surface the frequency ω decreases as the spinner constant α increases (figure 5a).
(ii) The low-frequency range
We examine the behaviour of ω(k) in the low-frequency limit for small values of |k| and evaluate the effective group velocity, corresponding to a quasi-static response of a homogenized elastic solid. The solutions to equation (2.15) may be expanded for small values of k1 and k2 as follows: 2.16 2.17
For the case when no spinners are present, these formulae reduce to
These expressions for and are consistent with the special case of eqn (3.3) in Colquitt et al. (2011).
In the supercritical regime m2≤α2, for low frequencies, only one dispersion surface exists representing a shear wave. For a monatomic, harmonic scalar lattice (with no spinners), the dispersion equation is given in Ayzenberg-Stepanenko & Slepyan (2008) as where the subscript s corresponds to the scalar case of out-of-plane shear. The low-frequency approximation for this dispersion equation is 2.18
It is interesting to compare equations (2.16)–(2.18) in the regime m2≤α2. Because of the presence of only one dispersion surface in this regime for the vector lattice of elasticity, this may be compared with an ‘equivalent’ scalar lattice in which the waves have the same dispersion properties as the shear waves in the vector lattice in the low-frequency regime. Comparing equations (2.16) and (2.18) one can deduce the relation between the stiffness parameters cs, c, the masses ms, m, and the spinner constant α, so that the shear waves in both types of lattices have the same dispersion properties. To achieve that, the required relation takes the form
(iii) Bi-atomic lattice of the vortex type
Now we show the computations for Bloch–Floquet waves for the bi-atomic triangular lattice described in §2a. The dispersion equation is given by relation (2.12), where and Additionally, the stiffness matrix can be expressed in the form (see Martinsson & Movchan 2003) with the 2×2 matrices and where ζ and ξ are given in equation (2.14).
For such a lattice with no spinners, the dispersion surfaces for Bloch–Floquet waves are shown in figure 6a (a section of the surfaces for k1=k2 is added in figure 6c). In contrast, for the vortex type of lattice, the corresponding results are given in figure 6b (with a section for k1=k2 in figure 6d).
These computations show that the vortex type of interaction within the lattice leads to formation of additional total band gaps, as shown in figure 6b. It is noted that the change in the behaviour of the lowest two dispersion surfaces through the introduction of spinners is small, whereas additional band gaps are introduced in the upper two dispersion surfaces.
The influence of the spinner constant α on the dispersion properties of Bloch–Floquet waves is illustrated in figure 7, where ω is represented as a function of α and kl=k1l=k2l. The two lower shear dispersion surfaces, represented in figure 7a, are defined for all real values of α, whereas the two upper pressure dispersion surfaces, depicted in figure 7b, disappear when critical regimes determined by the spinner constant α are reached. In particular, in the subcritical regime (α≤1), two pressure waves propagate while one disappears in the intercritical regime (1<α≤10) and no pressure waves propagate in the supercritical regime (α>10). Also, for the lower dispersion surfaces, at a sufficiently large α, the frequency ω decreases as the spinner constant α increases (figure 7a).
3. Continuous medium. Dynamic shielding
Here, we discuss the dynamic response and shielding properties of a continuous elastic solid whose equations of motion contain the vortex-type term.
This section is illustrative and is motivated by vectorial equation (2.3) in the long-wave approximation (|k|l≪1) for the monatomic triangular lattice which can be homogenized (see Gibson & Ashby 1988; Slepyan 2002). Here, we consider an elastic material with the Lamé constants λ=μ, mass density ρ and an additional term, representing the chiral properties of the medium, linked to a periodic system of spinners. The amplitude vector U of the time-harmonic displacement satisfies the full vector equations of motion in the plane (x,y) as follows: where μ is the shear modulus of the solid, Σ is the vorticity matrix defined in (2.13) and F is the body force. In the computations presented here, we normalize the mass density and the shear modulus of the chiral medium, so that ρ=1.0,μ=1.0, if not otherwise stated, and distinguish between (i) the subcritical case of |α|<1 and (ii) the supercritical case of |α|≥1, where α is the spinner constant from (2.13). Such a classification is linked to dispersion equation (2.15). In regime (I), there are two conical surfaces for small values of ω in the dispersion diagram (figure 3). On the other hand, in regime (II), there is just one acoustic dispersion surface, which has a conical shape for small values of ω.
(a) Point force in the chiral solid
The computations discussed in this section have been carried out using the COMSOL multi-physics finite-element software. To prevent reflection of waves from the artificial boundaries of the computational window, perfectly matched layer conditions for plane elasticity have been implemented (as in Brun et al. 2009).
Firstly, we show the results of the computational experiment produced for a continuous chiral solid loaded by a time-harmonic horizontal point force of the radian frequency ω=50. We do not relate these computations directly to the lattice model in the homogenization sense. Instead, we show the effect of chirality inducing shielding properties of the circular coating.
Figure 8 shows the change in the wave pattern for different values of the spinner constant α. In particular, figure 8a has the result of the computation for the non-chiral medium (α=0), while figure 8b includes the modified plot corresponding to a subcritical value of the spinner constant (α=0.6) demonstrating a vortex around the point force. Domination of pressure waves with a relatively large wavelength is clearly visible in the horizontal direction in figure 8a compared with figure 8b.
Figure 8c shows the critical case, with the spinner constant α=1.0, and the presence of spinners leads to a nearly isotropic map; the preferential directions associated with the point force are less pronounced.
Furthermore, figure 8d includes the displacement amplitude for another supercritical case with the larger value α=1.4, and the dominance of shear waves with the relatively small wavelength is clearly visible.
The computations suggest that the directional preference is substantially reduced in the presence of ‘spinners’ in the chiral medium in the supercritical regime of |α|≥1 and that the spinners induce shear polarization.
(b) An elastic inclusion with a chiral coating
Here, we consider an elastic plane containing a coated inclusion. The two-phase coating of the inclusion is assumed to have the chiral term in the governing equations, whereas the ambient medium around the coating and the inclusion itself are assumed to be non-chiral linear elastic isotropic materials. The coating consists of two semi-rings, with the spinner constants being of the same magnitude and opposite signs, as shown in figure 9. The coated inclusion is chosen to be symmetric about the straight line passing through the centre of the inclusion and the point where the concentrated force is applied.
An incident wave is generated by a point source placed at a finite distance from the inclusion. For the computational examples, the source term is provided by a point force or concentrated moment with a time-harmonic amplitude of radian frequency ω=50. In the computation, the elastic inclusion has the normalized Lamé constants λ=23 and μ=12 and the same normalized mass density ρ=1 as in the ambient medium. The coating is assumed to have the same elastic moduli λ=μ=1 and the same mass density as in the ambient elastic matrix.
For the isotropic medium with an inclusion without a coating (α=0), the results of the computation are shown in figure 10a. Comparing with figure 8a, one can clearly see the shaded region associated with the perturbation field produced by the elastic inclusion.
The next set of computations in figure 10b,c corresponds to the chiral two-phase coating placed around the inclusion and can be interpreted as a ‘cloak’ guiding shear waves around the inclusion. It is noted that, in these two diagrams where the point force is horizontal, the inclusion is placed in the region dominated by shear waves (characterized by relatively small wavelengths) It is also seen that the effect of cloaking is most pronounced in the supercritical regime where the spinner constant of the inclusion is chosen such that |α|≥1.
One of the important features of the computations shown in figure 10 is the directional preference of the ambient field produced by the point force. In particular, the aforementioned diagrams correspond to the case when the ambient field in the immediate neighbourhood of the inclusion is dominated by shear waves. A change of the orientation of the force with respect to the position vector of the inclusion results in a different pattern of the ambient field surrounding the inclusion. Figure 11 corresponds to the situation when the force is pointed towards (or opposite) the centre of the elastic-coated inclusion. As before, we show the results of the computations for the inclusion without the coating (α=0) in figure 11a, and supercritical chiral coating (α=1.5) in figure 11b. The direct comparison of these diagrams shows that the chiral coating reduces the shaded region associated with the inclusion. Furthermore, we observe that the ripples behind the inclusion have a size corresponding to that of the wavelength of shear waves, which suggests that the chiral coating acts as a shear polarizing cloak.
Finally, we consider the configuration when the ambient elastic field has no directional preference.
The numerical simulation in figure 12 shows the interaction of a radially symmetric shear wave, created by a time-harmonic concentrated moment of a radian frequency of ω=50, with an elastic inclusion; figure 12a presents the displacement amplitude in the plane with a disc-shaped inclusion without a coating (α=0); figure 12b shows the high level of suppression of the scattered fields when a coated inclusion (with a chiral coating with α=30) is introduced in the homogeneous medium. The graphs in figure 12 show the amplitude of the displacement measured along the main diagonal of the square domain connecting the upper left corner to the bottom right corner of the domain. The displacement behind the inclusion has been suppressed to a large extent. Needless to say, the applications of this type of design are limitless; for example, design of earthquake-resistant structures and drilling systems in geophysics.
The computation for the homogeneous medium without an inclusion is presented in figure 13a. We also note the results shown in figure 12a where the computation is presented for an uncoated inclusion. The shaded region can be reduced via introduction of the chiral coating. The supercritical case of α=2.0 is included in figure 13b, which shows a new diffraction pattern, with conical regions where waves have enlarged amplitudes. The two graphs in figure 13 show the amplitude of the total displacement along the main diagonal of the square computational region. They show that, along this line, the chiral coating has the capability to reconstruct the displacement pattern behind the inclusion in terms of amplitude and wavelength (typical of a shear wave).
4. Concluding remarks
For the first time, construction of an elastic metamaterial structure that possesses cloaking properties owing to in-built micro-rotations has been addressed for full vector, two-dimensional elasticity. The vorticity constants in the governing equations describing the chiral media are evaluated explicitly and, furthermore, dynamic responses of both discrete and continuous systems have been analysed in detail.
Novel properties of Bloch waves in discrete vortex-type elastic lattices have been identified with analytical findings complemented by numerical illustrations of dispersion surfaces.
In the continuum model, the chiral material has been used to design a composite cloak around an inclusion, so that the incident waves are guided around the inclusion. Such an interaction can lead to a substantial reduction of the shaded region behind the inclusion and to the observation that the coating can be interpreted as a polarizing cloak.
Furthermore, it was demonstrated that the parameters of the chiral coating can also be tuned so that the shaded region is enhanced and the amplitude of the displacement behind the inclusion becomes negligibly small.
We envisage a range of applications for the proposed model in the fields of geophysics and structural design of cloaks shielding defects, which interact with elastic waves.
This research was carried out while M. Brun was at the University of Liverpool as a Visiting Scientist. The financial support from the Research Centre in Mathematics and Modelling of the University of Liverpool and Italian Ministry of University and Research under PRIN 2008 ‘Complex materials and structural models in advanced problems of engineering’ is gratefully acknowledged.
↵1 Normalization has been used throughout the text, so that c=1, m=1 and the distance between the neighbouring masses is equal to unity. All other physical quantities have been normalized accordingly, and the physical units are not shown.
- Received March 16, 2012.
- Accepted May 1, 2012.
- This journal is © 2012 The Royal Society