Contents lists available at ScienceDirect

International Journal of Mechanical Sciences

journal homepage: www.elsevier.com/locate/ijmecsci

A semi-analytical approach towards plane wave analysis of local resonance metamaterials using a multiscale enriched continuum description

A. Sridhar, V.G. Kouznetsova*, M.G.D. Geers

Eindhoven University of Technology, Department of Mechanical Engineering, P.O. Box 513, Eindhoven 5600 MB, The Netherlands

CrossMark

A R T I C L E I N F O

A B S T R A C T

Keywords: Local resonance Acoustic metamaterials Enriched continuum Semi-analytical Multiscale Acoustic analysis

This work presents a novel multiscale semi-analytical technique for the acoustic plane wave analysis of (negative) dynamic mass density type local resonance metamaterials with complex micro-structural geometry. A two step solution strategy is adopted, in which the unit cell problem at the micro-scale is solved once numerically, whereas the macro-scale problem is solved using an analytical plane wave expansion. The macro-scale description uses an enriched continuum model described by a compact set of differential equations, in which the constitutive material parameters are obtained via homogenization of the discretized reduced order model of the unit cell. The approach presented here aims to simplify the analysis and characterization of the effective macro-scale acoustic dispersion properties and performance of local resonance metamaterials, with rich micro-dynamics resulting from complex metamaterial designs. First, the dispersion eigenvalue problem is obtained, which accurately captures the low frequency behavior including the local resonance bandgaps. Second, a modified transfer matrix method based on the enriched continuum is introduced for performing macro-scale acoustic transmission analyses on local resonance metamaterials. The results obtained at each step are illustrated using representative case studies and validated against direct numerical simulations. The methodology establishes the required scale bridging in multiscale modeling for dispersion and transmission analyses, enabling rapid design and prototyping of local resonance metamaterials.

© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Acoustic metamaterials can be used to engineer systems that are capable of advanced manipulation of elastic waves, such as band-stop filtering, redirection, channeling, multiplexing etc. which is impossible using ordinary materials [1-4]. This can naturally lead to many potential applications in various fields, for example medical technology, civil engineering, defense etc. The extraordinary properties of these materials are a result of either one or a combination of two distinct phenomena, namely local resonance and Bragg scattering. Bragg scattering is exhibited by periodic lattices in the high frequency regime where the propagating wavelength is of the same order as the lattice constant. Local resonance, on the other hand, is a low frequency/long wavelength phenomena and, in general, does not require periodicity. "Local resonance acoustic metamaterial" (LRAM) is therefore the term used here to distinguish the subclass of acoustic metamaterials based solely on local resonance. The present work is only concerned with the modeling and analysis of LRAMs restricted to linear elastic material behavior in the absence of damping.

* Corresponding author. E-mail address: v.g.kouznetsova@tue.nl (V.G. Kouznetsova).

Based on the primary medium of the wave propagation, LRAMs can be further classified as solid (e.g. [5]) or fluid/incompressible media (e.g. [6]) based. This paper is concerned only with the modeling of solid type LRAMs. A typical representative cell of such LRAMs is characterized by a relatively stiff matrix containing a softer and usually heavier inclusion. The micro-inertial effects resulting from the low frequency vibration modes of the inclusion induce the local resonance action. The complexity of the geometry of the inclusion plays an important role in the response of the LRAM. For instance multi-coaxial cylindrical inclusions have been proposed [7], which exhibit more pronounced micro-inertial effects due to the larger number of local resonance vibration modes, hence leading to more local resonance bandgaps. The symmetries of the inclusion also play a key role. A 'total ' (or omni-directional) bandgaps, where any wave polarization along any given direction is attenuated within a given frequency range is observed in micro-structures exhibiting mode multiplicity (or degenerate eigenmodes) resulting from a combined plane and 4-fold (w.r.t. 90° rotation) symmetry. For geometries with only plane symmetry, 'selective' (or directional) bandgaps, which only attenuate certain wave modes in a given frequency range can be observed. Furthermore, if the inclusion is not plane symmetric with respect to the direction of wave propagation, hybrid wave modes, which are a combination of compressive and shear wave modes can be

http://dx.doi.org/10.1016/j-ijmecsci.2017.08.027

Received 20 April 2017; Received in revised form 4 August 2017; Accepted 10 August 2017 Available online 13 August 2017

0020-7403/© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. (http://creativecommons.org/licenses/by/4.0/)

observed. In general, a more extensive micro-dynamic behavior results from an increased complexity in the LRAM design.

A plethora of approaches are available for modeling LRAMs. For steady state analysis of periodic structures, the Floquet-Bloch theory [8] gives the general solution that reduces the problem to the dispersion analysis of a single unit cell. The unit cell problem can then be further discretized and solved using appropriate solution methodologies. A popular technique highly suited to composites with radially symmetric inclusions is given by Multiple Scattering theory (MST) [9,10]. However, MST has not been successfully applied to more complex unit cell designs. The finite element (FE) method provides a robust approach for modeling arbitrary and complex unit cell designs. Such an approach, also known as the wave finite element (WFE) method, has been extensively employed in literature for analyzing periodic structures including acoustic metamaterials [11-14].

Since LRAMs operate in the low frequency, long wave regime, it is appropriate to exploit this fact towards approximating the general solution, leading to simpler methodologies. This is equivalent to ignoring nonlocal scattering effects resulting from reflections and refractions at the interfaces of the heterogeneities. This is the formal assumption made in many dynamic homogenization/effective medium theories [1521] which recover the classical balance of momentum, but where the local resonance phenomena manifests in terms of dynamic (frequency dependent) constitutive material parameters. The most striking feature of the effective parameters is that they can exhibit negative values, which indicates the region of existence of a bandgap. The specific macro-scale coupling of the local resonance phenomena is dictated by the vibration mode type of the inclusion [15]. A negative mass density is obtained in LRAMs exhibiting dipolar resonances e.g., [22]. Similarly, a negative bulk and shear modulus are obtained in LRAMs exhibiting monopolar and quadrupolar resonances, e.g. [23], though not further considered in this work. The homogenized models accurately capture the dispersion spectrum of LRAMs provided that the condition on the separation of scales is satisfied, i.e. the wavelength of the applied loading is much larger compared to the size of the inclusion. This can be ensured in the local resonance frequency regime by employing a sufficiently stiff matrix material compared to that of the inclusion.

The expression for the effective parameters of materials with radially symmetric inclusions has been obtained using analytical methods using the Coherent Potential Approximation (CPA) approach [15,16]. A generalization towards ellipsoidal inclusions has also been presented in [17,24], thereby extending the method of Eshelby [25] to the dynamic case. Analytical models for unit cells composed discrete elements (e.g. trusses) have also been derived [19,26,27]. For arbitrary complex micro-structures, it is necessary to adopt a computational homogeniza-tion approach. A FE based multiscale methodology in the framework of an extended computational homogenization theory was presented in [20,21].

The present work builds upon the computational homogenization framework introduced in [21]. A first-order multiscale analysis is combined with model order reduction techniques to obtain an enriched continuum model, i.e. a compact set of partial differential equations governing the macro-scale behavior of LRAMs. The reduced order basis is constructed through the superposition of the quasi-static and the micro-inertial contribution, where the latter is represented by a set of local resonance eigenmodes of the inclusion. In the homogenization process, the generalized amplitudes associated to the local resonance modes emerge as additional kinematic field quantities, enriching the macro-scale continuum with micro-inertia effects in a micromorphic sense as initially defined by Eringen [28]. An equivalent approach for modeling LRAM was derived using asymptotic homogenization theory in [18], but was not further elaborated and demonstrated as a computational technique.

In this work, the enriched continuum of [21] is exploited to develop an ultra-fast semi-analytical technique method for performing dispersion and transmission analysis of LRAMs with arbitrary complex microstructure geometries that retains the accuracy of WFE methods at a frac-

tion of the computational cost. A modified transfer matrix method is developed based on the enriched continuum that can be used to derive closed form solutions for wave transmission problems involving LRAMs at normal incidence. Furthermore, the dispersion characteristics of the enriched continuum and their connection to various symmetries of the inclusion are discussed in detail. In order to ensure the accuracy of the method in the frequency range of the analysis, a procedure for verifying the homogenizability of a given LRAM is introduced. Although the enriched continuum predicts both negative dynamic mass and elastic modulus effects, the latter exhibits a significant coupling only at frequencies close to and beyond the applicability (homogenizability) limit of the present approach. Since such effects cannot be robustly modeled, they will not be considered in the present paper.

The structure of the paper is as follows. The relevant details of the enriched multiscale methodology are briefly recapitulated in Section 2. In Section 3, a plane wave transform is applied to obtain the dispersion eigenvalue problem of the enriched continuum. The influence of the inclusion symmetry on the dispersion characteristics is highlighted. Based on the obtained dispersion spectrum, a procedure for checking the applicability (homogenizability) of the problem in the frequency range of analysis is elaborated that provides a reasonable estimation of its validity. In Section 4, a general plane wave expansion is applied to derive a modified transfer matrix method for performing macro-scale transmission analyses of LRAMs at normal wave incidence. The results obtained in each of the sections are illustrated with numerical case studies and validated against direct numerical simulations (DNS). The conclusions are presented in Section 5.

The following notation is used throughout the paper to represent different quantities and operations. Unless otherwise stated, scalars, vectors, second, third and fourth-order Cartesian tensors are generally denoted by a (or A), a, A, A(3) and A(4) respectively; ndim denotes the number of dimensions of the problem. A right italic subscript is used to index the components of vectorial and tensorial quantities. The Einstein summation convention is used for all vector and tensor related operations represented in index notation. The standard operations are denoted as follows for a given basis ep, p = 1, . .ndim, dyadic product: a®b = apbqep ® e, dot product: A • b = Apqbqep and double contraction: A : B = ApqBqp. Matrices of any type of quantity are in general denoted by (•) except for a column matrix, which is denoted by (;) . A left superscript is used to index quantities belonging to a group and to denote sub-matrices of a matrix for instance for a and by mn a and respectively. The transpose of a second order tensor is defined as follows: for

A = Apq ep (

eq, AT = Agpep ® eq. The transpose operation also simulta-

neously yields the transpose of a matrix when applied to one. The first and second time derivatives are denoted by (.) and (•) respectively. The zero vector is denoted as 0.

2. An enriched homogenized continuum of local resonance metamaterials

This section summarizes the essential features of the enriched continuum model introduced in [21]. The framework is based on the extension of the computational homogenization approach [29] to the transient regime [20]. The classical (quasi-static) homogenization framework relies on the assumption of a vanishingly small micro-structure in comparison to the macroscopic wavelength. No micro-inertial effects are recovered in this limit. The effective mass density obtained in this case is a constant and is merely equal to the volume average of the microstructure densities, unlike the frequency dependent quantity observed in metamaterials.

The key aspect rendering the classical computational homogeniza-tion method applicable to LRAMs is the introduction of a relaxed scale separation principle. Making use of the typical local confinement of the resonators in LRAM, the long wave approximation is still assumed for the matrix (host medium), while for the heterogeneities (resonators), full dynamical behavior is considered. Note, that the long wavelength

approximation poses a restriction on the properties of the matrix which has to be relatively stiff in order to ensure that it behaves quasi-statically in the local resonance frequency regime.

With this relaxed separation of scales, the multiscale formulation for LRAM is established. The full balance of linear momentum equations are considered at both the micro and macro-scales. The domain of the micro-scale problem associated to each point on the macro-scale domain, termed the Representative Volume Element (RVE), is selected such that it is statistically representative, i.e. captures all the representative micro-mechanical and local resonance effects at that scale. For a periodic system, the RVE is simply given by the periodic unit. The homogenization method can be applied to a random distribution of inclusions provided that there is sufficient spacing between the inclusions (i.e. excluding the interaction between the resonators).

The Hill-Mandel condition [30], which establishes the energy consistency between the two scales, is generalized by taking into account the contribution of the momentum into the average energy density at both scales. As such, the resulting computational homogenization framework is valid for general non-linear problems. Restricted to linear elasticity, a compact set of closed form macro-scale continuum equations can be obtained. To this end, exploiting linear superposition, the solution to the micro-scale problem is expressed as the sum of the quasi-static response, which recovers the classical homogenized model, and the micro-inertial contribution spanned by a reduced set of mass normalized vibration eigenmodes s<°, of the inclusion (shown in continuum form) as follows,

"m (*m) = " + §(3) (*m) : V" + £ "$(xm)% , (1)

where um and xm represent the micro-scale displacements and position vector respectively, ÏÏ the macro-scale displacement and stj the generalized (modal) amplitude of the s^1 local resonance eigenmode of the inclusion indexed by the set Q. The size of Q denoted by N(Q) gives the number of modal degrees of freedom. The third order tensor S(3) (3m ) represents the quasi-static response of the RVE to applied macro-scale strain under periodic boundary conditions. The eigenmodes are obtained via an eigenvalue analysis of a single inclusion unit under fixed displacement boundary condition. For LRAMs, these eigenmodes naturally represent the local resonance vibration modes of the system. Only one inclusion needs to be considered here even for random distributions since all relevant micro-inertial properties can be obtained from it. The generalized degrees of freedom associated to these eigenmodes then emerge at the macro-scale as additional internal field variables, thus resulting in an enriched continuum description. FE is used to discretize and solve for the quasi-static response and the eigenvalue problem.

Further details of the derivation of this framework can be found in [21]. In the following, only the final equations describing the homogenized enriched continuum are given.

The macro balance of momentum

V • a1 - p = 0 ) (2)

The reduced micro balance of momentum (to be solved at the macro-scale)

s®?es t + t = -SJ •ÏÏ-SH :VH,se Q) (3)

The homogenized constitutive relations

g = C (4) : Vo + — £ SH stt, (4a)

— seQ

P = PU + -1 ^ sjst. (4b)

— seQ

Here, a represents the macro-scale Cauchy stress tensor and p the

— (4) _

macro-scale momentum density vector. The terms C and p are, respectively, the effective linear (static) elastic stiffness tensor and the

effective mass density of the RVE (the over-bar (5) is used to distinguish an effective material property from its respective counterpart of a homogeneous material). These are the classical quantities obtained under the quasi-static approximation of the micro-structure and are not to be confused with their dynamic counterparts [15]. The term swres represents the eigenfrequency (in radians per second), of the sh local resonance eigenmode of the inclusion. The formulation of depends on the minimum number of eigenmodes, required to sufficiently accurately capture the dispersive behavior of the system in the desired frequency regime. A more precise mode selection criterion is discussed in Section 3. The set Q is indexed in the order of increasing eigenfrequencies, i.e. for r, s eQ, r < s implies rrnles < scores. The coupling between the macro and reduced micro-scale balance equations is represented by the vector sj and the tensor SH. The vector sj describes the coupling of dipolar local resonance modes, whereas SH gives the coupling of monopolar, torsional and quadrupolar modes. Finally, in Eq. (4), V represents the volume of — (4) _ " —

the RVE. The coefficients C , p, sj, SH and s®res (s e Q) constitute the set of effective material parameters characterizing an arbitrary LRAM RVE. They are obtained through discretization and model reduction (See [21] for details). Specifically, the parameters sj and SH are obtained by projecting the computed inertial force per unit amplitude of the local resonance eigenmodes onto the rigid body and the static macro-strain deformation mode respectively. The explicit expressions for these effective parameters in continuum form are given as follows,

SJ = [ P$m) sMn ) dV, (5a)

SH = Í sk*m) • (p(v ) §(3) (2m ) )d V, (5b)

where p( xm ) represents the mass density of the heterogeneous microstructure. Subsequently, the system (2)-(4) can be solved at the macro-scale.

In the following sections, the analysis will be restricted to dipolar local resonances, i.e. the coefficient SH will be disregarded in the sequel. The justification for this is as follows. Since the primary concern of the present paper is to develop an accurate methodology comparable to standard numerical methods, the analysis is restricted to the deep sub-wavelength regimes, where the applicability of the method (i.e. ho-mogenizability of the problem) is guaranteed. The details of the homog-enizability criterion are elaborated in Section 3.3, where it is required for the effective stiffness of the matrix structure (both shear and com-pressive) to at least 100 times higher than that of the inclusion. In this case, the boundary of the inclusion is effectively rigid, thereby suppressing the action of monopolar and quadrupolar eigenmodes, which necessarily require a deformable boundary [15]. This fact has also been evidenced in [18] where a homogenized model has been derived via asymptotic homogenization. Indeed, within the homogenizability limit of the proposed method, the bandgaps due to negative stiffness effects have negligible bandwidths and appear as flat branches (see for e.g. the branch corresponding to mode 7 in Fig. 2a representing a monopolar resonance). Since only dipolar resonances are dominant in the regime considered here, it is therefore justified to neglect the term SH in further derivations for the sake of simplification.

3. Dispersion analysis

In this section, a plane wave analysis is carried out on an infinite enriched continuum, i.e. without taking into account macro-scale boundary conditions. First, a plane wave (Fourier) transform is applied on Eqs. (2)-(4) to obtain the dispersion eigenvalue problem valid for an arbitrary LRAM RVE design. The relation between the emergent dispersion characteristics and the properties of the dynamic mass density tensor and the geometric symmetries of the inclusion is then discussed. A procedure for checking the homogenizability of the approach is then

■ Lead

□ Soft Rubber

□ Hard Rubber

□ Epoxy

Fig. 1. The unit cell designs used for the case studies (a) UC1, (b) UC2 and (c) UC3.

established via dispersion analysis. Three example periodic unit cell designs are introduced as case studies to illustrate the influence of the inclusion geometry on the general dispersion characteristics of LRAMs and to validate the methodology against standard Bloch analysis.

3.1. Theoretical development

The plane wave transform of any continuous field variable can be expressed as,

where the (J) represents the transformed variable, k denotes the wave vector, m the frequency and i the imaginary unit. The wave vector can be represented by its magnitude and direction as k = keg, where k = is the wave number, X the corresponding wavelength and eg a unit vector in the direction of wave propagation. The above expression provides the plane wave transform in an infinite medium. Applying the transformation (6) to Eqs. (2)-(4) while disregarding SH as discussed before (see the last paragraph in Section 2) yields,

ik ■ a + imp = 0 ,

S®?es Sf> - °°2 Sf> = ®2 SJ ■ s e Ö

I-(4) - \ p

il ( C) ■ k) ■u,

p = — irnt pu +

(7a) (7b) (7c)

Eliminating a and p by substituting Eq. (7c) and (7d) into equation (7a) and combining it with Eq. (7b) results in the following parameterized eigenvalue problem in m,

k2 Ce Q

o>2 —res

r - r -• ]

■u 0

• ®res is a Nq x Nq diagonal matrix of the eigen-frequencies, tf and j are column matrices of size Nq containing the

modal amplitudes s tj and the coupling vectors sj, respectively; Q and O are column matrices of size Nq with all entries equal to the zero vector 0 and scalar 0, respectively. The above eigenvalue problem is the so-called k - m form, in which m is the eigenvalue and k and ee are parameters. The problem consists of ndjm + Nq variables, which gives the number of dispersion branches predicted by the model. The alternative m - k form in which k is the eigenvalue and m and ee are the parameters can be obtained by eliminating tf from the system of Eqs. (8) resulting in

[k2 Ce — a? where

p(m) ) ■u =

p(ffl) = p\ + ^

where SJ = ^ sj®sj is a measure of the translational inertia associated to each local resonant eigenmode, termed modal mass densities and ~p(ai) is the effective dynamic mass density tensor which is a well known quantity in the literature [18,22,31,32,33]. In most works, however, the expression for the dynamic mass density is derived for a particular unit cell design often after some simplifications and approximations. Here, on the contrary, Eq. (10) holds for arbitrary geometries, for which it is obtained numerically.

The solution to Eq. (9) at a given m yields Udim eigen wave numbers kp(m), p = 1. .ndjm and corresponding eigen wave vectors denoted by vp(m). These vectors are also termed polarization modes. Let vp(m) be normalized with respect to Ce, i.e.

op ■ Ce ■ i

In the present context, only positive real or imaginary eigenvalue solutions are considered to avoid ambiguity in the direction of wave propagation or decay respectively. If vp(m) is parallel to ee (i.e. vp x ee = 0) for a given p , the corresponding wave mode is purely compressive. On the other hand if vp(m) is orthogonal to ee (i.e. ee • vp = 0), the resulting mode is purely shear. All other combinations represent hybrid modes.

For an arbitrary RVE geometry and wave direction, the dispersion problem (8) and (9) can be solved numerically which, owing to the reduced nature of the enriched continuum problem, is computationally much cheaper in comparison to standard Bloch analysis techniques, especially in the 3D case. It is arguably cheaper even compared to reduced order Bloch analysis techniques [34,35] since the spectrum of the homogenized problem is much smaller than the Bloch spectrum leading to a more compact dispersion equation.

If accounts for all dipolar modes, the following relation holds due to the mass orthogonality of eigenmodes

«nc P!>

where ^inc is the static mass fraction of the inclusion defined as the ratio of the mass of the inclusion (including the central mass and the compliant support) and the total mass of the RVE. The above relation can be used to formulate a mode selection criterion for setting up . Projecting Eq. (12) along some arbitrary set of orthonormal basis vectors ep, p = 1. .ndjm gives a set of scalar equations,

y ßjp = ßinc -see

seö r es

where, s^Jp = (p)"1 ep • SJ • ep, s e Q is the modal mass fraction of the local resonant eigenmode along ep for p = 1. .ndim. In practice, it is sufficient if the above expression is satisfied in the approximation as the local resonance modes in general will never fully be able to capture all the mass of the inclusion. The tolerance can vary based on the design requirements but a good measure would be 5% of the inclusion mass fraction. Thus, Q should be the smallest set of low frequency eigenmodes that satisfies Eq. (13) within the given tolerance.

ik^x—iat

— CO

s„ ,2

Fig. 2. The dispersion spectra of the considered unit cells, (a) UC1, (b) UC2 and (c) UC3 computed using Bloch analysis (black dashed lines) and the presented semi-analytical (SA) approach (colored lines). The mode shapes of some of the local resonant eigenmodes and their associated dispersion branches are also shown. The topmost colorbar represents the cosine of vp, p = 1, 2, with respect to e0, red indicating a pure compression wave and blue a pure shear wave, intermediate colors represent hybrid wave modes. The colorbar below indicates the norm of the displacements of the local resonance eigenmodes. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3.2. Relation between the dispersion characteristics, the dynamic mass density tensor and geometrical symmetries of the inclusion

The dispersion characteristics of a LRAM are primarily determined by the properties of ~p(a>) . It is anisotropic in general and its components can take all values from [—<x, at different frequencies. A bandgap

exists at m when at least one of the eigenvalues of ~p(a>) is negative. From Eq. (10), it can be concluded that a bandgap is always initiated at sm, s eQ , where the determinant of ~p jumps from » to —». For frequencies close to s m, the dispersion is strongly determined by the properties of the corresponding sj vector due to the large magnitude of its coefficient in Eq. (10). Around that frequency, the polarization modes vp, p = 1, . .ndim obtained from the dispersion problem (9) become

oriented with respect to sj, either along the vector or orthogonal to it. If sj is not parallel or orthogonal to ee, hybrid wave solutions are obtained. Note, that hybrid wave modes can also result from the anisotropy of the effective static stiffness tensor. However, only the hybrid modes resulting from local resonances are considered here.

A "total" bandgap is formed in the frequency range where all the eigenvalues of ~p are negative (i.e. it is negative definite). Within this range, all waves regardless of the direction will be attenuated. When ~p is only negative semi-definite, a "selective" bandgap is formed that attenuates some of the wave modes. Two bandgaps will overlap only if a pair of eigenfrequencies s m, r m, s, r e Q, s ^r are sufficiently close enough and sj is orthogonal to rj. In the special case of a system exhibiting mode multiplicity, i.e. sm = rm and ||sj ||= \\rj ||, the bandgaps will overlap perfectly. Such a situation is guaranteed for an inclusion domain with isotropic constituents that possesses the following symmetries; plane symmetry with respect to 2 mutually orthogonal axes, and a 4-fold (discrete 90°) rotational symmetry about the normal to these axes (exhibited for e.g. by a square, cross, cylinder, hexagonal prism etc.). In this case, all rj vectors will be aligned along one of the symmetry axis. For a 3D inclusion, a mode multiplicity of 3 is observed provided the inclusion has plane symmetry about 3 mutually orthogonal axes and 4-fold rotational symmetry about these axes.

For such a highly symmetric inclusion domain, the set can be re-parameterized by gathering all the sj vectors with duplicate eigenfrequencies and labeling them using the same number. Let Qs represent

such a set. Here N (Q s ) :

N(Q). Let e , p = 1, ..ndim be the unit vec-

tors representing the symmetry axis. For every s e Qs, a set of vectors

Is 1 p = \[7JpVep, p = 1) ..ndim} is obtained where sn3 = (pV)-1 ||sjp||2, is the modal mass fraction associated to that degenerate eigenmode group. With this, the expression for the effective mass density tensor (10) for plane + 4-fold symmetry case can be re-written by gathering all the terms under the duplicate eigenfrequencies as follows,

p(m) = p\ + 2

=?(i + 2

■ (B2 V

seQs res

where the fact

'"dim £

p=1 P '

>ep = I has been used. Thus, the mass density

tensor becomes isotropic and diagonal in this case. Such a tensor is always either positive or negative definite at any given frequency. Therefore, an inclusion with plane and 4-fold symmetry only exhibits total bandgaps and no hybrid wave solutions.

For inclusions with only plane symmetry with respect to ep, p = 1 , . ."dim (but not the 4-fold symmetry), for e.g. an ellipsoid, rectangle etc., degenerate eigenmodes are no longer guaranteed, however the sj vectors will still align along ep. Now the set Q is divided into n^dim subsets pQ' defined as pQ' = js e Q| ep ■ SJ • ep ^ 0} for p = 1..ndlm. Using Eq. (10), the mass density tensor can now be re-written for the plane symmetry case as follows,

P(b) = Pi

i + 2 2

V P=1 sep Q.

( 2 u res

s ßj ep ® e p

where s j = (pV) 1 ||sj||2. Thus the dynamic mass density tensor is now orthotropic, or diagonal with respect to the axis defined by ep,

p = 1, . .ndim. If the wave propagates along one of the symmetry axis, i.e. eg = ep for some p = 1, . .ndim, then no hybrid wave mode solutions will be obtained at all frequencies due to local resonance. Furthermore, if Cg can be diagonalized with respect to ep, which is the case if the effective stiffness tensor is either isotropic or orthotropic with the orthotropy axis aligned with the symmetry axis of the inclusion, Eq. (9) can be fully decoupled giving ndim independent scalar equations,

k2pСврр - (в2 ppp(œ) = 0 , (no summation on p),

where,

Pvv(a>) = P\

1 + jQ

ч s&Q'

Here, Cgpp = ep • Cg • ep and ppp(m) = ep • p(m) • ep. The solution for the wave number derived from equations (16) and (17) is given as

kp(m) ■■

с0вр

(i + Z тт^^)'

\ S&Q> S®ies - J

where,

represents the effective wave speed in the quasi-static limit of the system. In accordance with the normalization of vp with respect to Cg given by Eq. (11) , the expressions for the wave polarization modes are

Hence, in the plane symmetric case, one compression mode (with vp || ee for a given value of p ) and two shear (one in the 2D case) modes (with vq 1 eg for p i=- q) are always observed without formation of hybrid wave modes, for wave propagation along the symmetry axis.

Since the combined plane and 4-fold symmetry is a special case of plane symmetry, the scalar dispersion relations (18) also hold in this

3.3. Homogenizability limit

A criterion on the applicability limit of the developed semi-analytical analysis can be obtained heuristically. The relaxed scale separation principle (stated above in Section 2) no longer applies when the wavelength of the macroscopic wave in the matrix approaches the size of the heterogeneities. At that limit, higher order scattering effects start to play a role which is not accounted for by the present homogenization theory. A safe estimate for the scale separation limit is when the matrix wavelength is at least 10 times the relevant microstructural dimension, i.e.

kp < 0. 2 j, where € is the relevant dimension for all p ■■

1, . ,nA

For the local resonance frequencies to lie within this limit, the effective stiffness of the matrix structure (both shear and compressive) has to be at least 100 times higher than the effective stiffness of the inclusion. The high stiffness of the matrix structure also implies that the local resonances will be internally contained, i.e. preventing any interaction between neighboring inclusions. This also imposes a constraint on the maximum volume fraction of the inclusions, which is largely determined by the material properties of the matrix and the desired frequency limit of analysis. A stiffer matrix material can internally contain the local resonances at higher volume fractions compared to a compliant matrix. Similar arguments apply to the cases of random inclusion distributions where the inclusions can be arbitrarily close to one another, which can ostensibly result in some interaction. Hence it is assumed here that, in the case of random inclusion distributions, there is still sufficient distance between any two neighboring inclusions.

The approximate frequency limit of the homogenization method can be estimated by solving the dispersion Eq. (8) at k = 0 . 2 p along a given

Table 1

The geometric and material parameters of the considered unit cell designs.

(a) Geometric parameters of the unit cell Dn Diameter of lead inclusion

Dout Outer diameter of rubber coating

Width of hard rubber insert Length of the unit cell = f2 ) Volume of unit cell

(b) Istropic linear elastic material parameters

Material

p [kg/m3 ] E [MPa]

Epoxy 1180

Lead 11600

Soft rubber 1300

Hard rubber 1300

3.6 x 103 4.082 x 103 0.1175 11.75

10 mm 15 mm 2 mm 20 mm

400 x 103 mm 3

0.368 0.37 0.469 0.469

wave direction. The smallest eigenfrequency for which the corresponding group velocity is non-negligible, indicating a matrix dominant wave, gives a conservative estimation of the homogenizability limit of the method. Local resonance modes occurring beyond this limit might exhibit significant nonlocal interactions between neighboring inclusions, at which point the present homogenization method will become inaccurate and ultimately fail.

(19) 3.4. Numerical case study and validation

In order to validate the proposed semi-analytical methodology and to illustrate the results presented thus far, the dispersion properties of three 2D LRAM RVE designs shown in Fig. 1 are studied. A periodic micro-structure is assumed, hence the RVEs considered here represent the periodic unit cell. The periodicity is assumed for the sake of simplicity of analysis since the distribution is captured by a single inclusion. It also allows Bloch analysis [14] to be performed in order to compute the reference solution. However, the results obtained for a periodic unit can to a large extent be applied to a random distribution of the inclusions as well, provided that the average volume of the matrix material over the macroscopic domain remains the same and that the inclusions are well spaced between each other, i.e. the localized resonances are not influenced. This is because the effective dynamic mass density tensor, which solely determines the bandgap properties (i.e. its position and size) resulting from local resonance, is unaffected by the distribution of the inclusions. The only difference between a periodic and a random distribution is the degree of anisotropy of the resulting effective static stiffness. Thus, only the phase speeds of the propagating waves along oblique angles with respect to the lattice vectors will be different in this case.

The considered unit cell designs are based on that of Liu et al. [5], each consisting of a square epoxy matrix with an embedded rubber coated cylindrical lead inclusion. The difference between the designs is in the configuration of the rubber coating such that each unit cell exhibits a different degree of symmetry. The first design in Fig. 1a is both plane symmetric with respect 0j and 02 and 4-fold symmetric about 03 . The second design in Fig. 1b is only plane symmetric with respect to the 0j - 03 and 02 - 03 planes, in which a patch of hard rubber replaces the soft rubber as shown in the figure. The third design in Fig. 1c is obtained by rotating the inclusion in the second design 22.5° counterclockwise with respect to 03 . These three examples will exemplify the influence of the dynamic mass anisotropy on the formation of total bandgaps, selective bandgaps and hybrid wave modes. The designs are henceforth labeled as UC1, UC2 and UC3, respectively. Plane strain kinematics is assumed in the 0j - 02 plane. The values of the geometric and material parameters of the three designs are given in Table 1a and Table 1b, respectively1. A finite element discretization is used to obtain the numerical models of the unit cells. Four-node quadrilateral finite elements were used with a maximum mesh size restricted to 0.57 mm in the matrix and 0.2 mm in the rubber coating. The models comprise on average

1 The material properties of the two rubber materials used here are hypothetical and do not target a particular rubber material.

Table 2

Homogenized enriched continuum material properties of the considered unit cells.

(a) Effective static material properties. Only non-zero components of C are shown. Design UC1 & UC2 & UC3

[GPa] Cirn = C2222 = 1. 955

C1122 = C2211 = 0 . 546

C1212 = C2121 = C1221

p [kgm-3] 3201.6

(b) Enriched continuum effective properties associated to the local resonance eigenmodes. Here e' and el° rotated by 22.5° with respect to e1 and e2, respectively Design UC1 UC2 & UC3

sœ [Hz]

r= 1 2

355.7 0.724

1239 0.0468

r = 2 3

355.7 0.724

1239 0.0468

•j_= i C1 = 1

508 0.75

' ßJ Pe1

/"J I' ) 6

1615 0.0157

•j_=

C j = 3

936 0.704

ßjP <¡2 ßjp ¡2) 4

1263 0.0675

2053 0.0013

8000 elements and 15,000 degrees of freedom. Convergence of the results with respect to the discretization size has been verified.

The homogenized enriched continuum parameters are extracted from the finite element models by applying the method described in [21]. The values of the computed static effective parameters are given in Table 2a. Due to the high compliance of the rubber coating in comparison to the epoxy, the overall effective stiffness of each unit cell only differs up to 1%, and hence this difference is neglected in the Table. The effective mass densities of all the designs are identical. The first 15 local resonance eigenmodes are determined and the corresponding enriched material parameters are computed for each case. The mode shapes of some of the modes are displayed in Fig. 2 for reference. Among these, only modes 2, 3, 5 and 6 of UC1 and modes 1, 3, 4, 6 and 9 of UC2 and UC3 unit cells are dipolar and possess a sufficient modal mass fraction along one of the considered directions. Modes 1 and 4 of UC1 represent the torsional resonances, involving the central inclusion and the coating, respectively. Mode 7 of UC1 is of monopolar type as evidenced from the fact that it is plane symmetric. The dipolar modes 2 and 3 in UC1 and 1 and 3 in UC2 and UC3 represent the vibration of the central core whereas modes 5 and 6 in UC1 and 4, 6 and 9 in UC2 and UC3 represent localized vibration in the coating without any involvement of the central inclusion. The properties of the dipolar modes are given in Table 2b. The sum of the modal mass fractions of the dipolar modes is almost 98% of the total mass fraction of the inclusion ( ninc = 0.787) along any given direction, hence forming a sufficient basis that satisfies the mode selection criterion given by Eq. (13). However, for the sake of validation, all 15 modes were considered in the analysis. Furthermore, similar to the static properties, the effective properties describing local resonance modes are identical for the UC2 and UC3 modes, with the only difference being the orientation of the coupling vector sj of UC3, which would be rotated by 22.5° counter-clockwise with respect to the corresponding ones of UC2.

The dispersion spectra for the considered unit cells are computed using Eq. (8) and validated against the spectra computed using Bloch analysis [14] as shown in Fig. 2. An excellent match is observed between the spectra computed using the two methods for all designs within the displayed frequency range. The group velocities observed at approximately 0.2 times the distance from r to X for the three unit cells are all close to zero, hence confirming again that within the considered frequency range, the analysis lies well within the homogenizability limit discussed in Section 3.3.

Only total bandgaps are observed in the case of UC1 as expected and mostly selective bandgaps are observed for UC2 and UC3. Two narrow total bandgaps are seen for UC2 and UC3 in a narrow frequency range above the eigenfrequncy of the 3rd and 6th eigenmode at 936 and 1615 Hz respectively as shown in Fig. 2b and 2c. The formation of hybrid wave modes is clearly visible in UC2 and UC3 along r-M and X-M, when the wave direction is not aligned along the respective symmetry axis of the inclusion. The hybrid wave effects are localized around the

first local resonance frequency. The loss of hybrid wave mode effects when the wave propagation is perfectly oriented along the symmetry axis can also be observed. Furthermore, the bandgaps of UC2 and UC3 appear at exactly the same frequency locations, which is a consequence of the fact that both designs possess the same set of eigenfreqencies and modal mass fractions. Several flat branches are also observed in both spectra corresponding to the modes with zero modal densities (i.e. tor-sional, monopolar modes etc.). Hence it is confirmed that these modes indeed do not influence the dispersive properties of the system and the modes indicated by the proposed selection criterion are adequate.

Finally, it should be emphasized that the total computational cost of the reduced dispersion problem including the offline cost of computing the enriched material parameters which involves meshing, matrix assembly and model reduction for a given unit cell design, is significantly smaller compared to Bloch analysis on the full FE model (even with all 15 modes included). The Bloch and the homogenized analysis were both executed using a MATLAB script on a desktop computer with an Intel i7-3630 QM core and 6GB memory. The Bloch analysis took on an average 230 s for each unit cell, while the homogenized analysis took a total time (including offline and online costs) of about 4 s, indicating a tremendous gain ( ~ 50 times in this case) in computational speed.

To summarize this section, the various local resonance effects resulting from complex micro-structure designs such as total and selective bandgaps and hybrid wave modes have been illustrated. The semi-analytical model accurately captures all effects in the given frequency regime, while being computationally much more efficient and faster compared to a full Bloch analysis.

4. Transmission analysis

In this section, a transmission analysis framework is derived for wave propagation in an enriched media at normal incidence. It is a generalization of the standard transfer matrix method [36] used in the analysis of plane wave propagation through a layered arrangement of several dissimilar materials, which can in general be distinct locally resonant metamaterials, described by enriched effective continuum. The technique can be extended to oblique wave incidence, but this is beyond the scope of the present paper. The framework is applied to analyze the steady state responses of systems constructed using PlySym and UC3 unit cell designs including the influence of the finite size of the structure and the boundary conditions.

4.1. Theoretical development

The general problem can be described as a serial connection of m layered (enriched) media with the first and the last layer being semi-infinite as shown in Fig. 3. Since only normal wave incidence is considered, the wave direction vector ee defines the 1D macro-scale coordinate axis. Let x represent the corresponding spatial coordinate. The general plane

forward backward wave wave

med 1. -► med 2.

med r.

r+l{ r+1 '

med r+1.

med m.

Fig. 3. The general macroscopic acoustic boundary value problem for normal wave incidence.

wave solution at x corresponding to medium r can be expressed as a sum of the forward and backward propagating components represented by % and % respectively. Hence,

u(x) = rUf(x) + rub(x).

The bounding semi-infinite media will posses only a forward wave or a backward wave depending whether its boundary is on the left or the right side, respectively. Each component of the total displacement is composed of the individual wave modes

actuator

homogeneous medium

Fig. 4. The acoustic transmission problem.

X = ni

f) = % %% e ' P=1

% (x) = 2 v^m

where, rf and r£bp, p = 1 , . .ndjm are the wave mode amplitudes of the forward and backward waves, respectively. These can be found by making use of the normalization condition (11), giving

f = ' Cg •r vF • %(x)e~'kx,

lb, = ' C • •'k( x)e'"k?x.

Next, the traction-displacement relation needs to be setup. Let rt+(x) denote the macro-scale traction vector in the r^1 medium with respect to eg (a " - " superscript is used to indicate traction with respect to -ee). This can be determined from the effective homogenized constitutive relation (4a), disregarding the last term as discussed before (see the last paragraph in Section 2), i.e.

rf (x) = ee • ra(x) = eg • "(C^ : V3(x).

The displacement gradient is expressed in terms of the plane wave solution by making use of Eqs. (21)-(23),

frequency dependent unlike the impedance of ordinary homogeneous elastic materials.

Another relation that is needed for solving the boundary problem is the transfer operation, which, given a quantity known at x, returns its value at another point x + Ax within the corresponding medium. The expression for components r(i)f( x + Ax) and r(J)b (x + Ax) in terms of

'(i)f(x) and '(i)b(x) respectively is obtained by applying Eqs. (22) and (23f), b

'(î)f(x + Ax) = £ ® 'Cg • ^ <AAx) • '(î)f(x) = <T(Ax) • f).

"dim __) _1

'(î)b(x + Ax) = £ (r«p ® rCg • rv„e-'rk"Ax) • r«b(x) = rT (Ax) • '(î)b(x).

where,

<T(Ax) = £ 'tp®'C • p =1

is the effective transfer operator.

Finally, the continuity conditions at the media interfaces are established. Let rx give the coordinate of the r^1 and (r + 1)'h interface. Applying Eqs. (21) and (26), the traction and displacement continuity conditions between the r and (r + 1)'h medium r x can respectively be written as,

^ „ dim , __^ _

VÜ(x) = fktY, ( ëg ®rUp® rCg •rDp • Vf(x) - U®rUp®rCg •r vp

% (x) f,

"dim __^ __„ )

'û\Ya T Cg • r°p • Vf(x) - eg ® r5p ®r Cg • rvp • ^(x) f,

where rKp(œ) = œ-1 rkp(œ) has been introduced. Substituting Eq. (25) in (24) gives,

r t+(x) = iœrZ(œ) • % - iœrZ(œ) • % , where,

rZ(œ) = £ ' Cg • (rKp(œ)r Df(œ) ® r Cg • ^ (œ)) ,

is the effective impedance of the considered metamaterial medium, defined as the constitutive parameter relating the macro-scale traction to the velocity ( iofiif or io/uhf at a given interface. Note from Eq. (26) that a positive sign is assigned to the impedance with respect to the forward wave component of the total velocity and a negative sign is assigned to the impedance with respect to the backward wave component. It is

r1+ ( fx) + r+1~t~ ( rx) = tt

iœrZ(œ) • CitfCx) - rÜbfrx)) - iœr+1Z(œ) • C+fx) - r+1fÜbffx)) = 0

Ü(rx) -

r Ü ff fx) + % ( fx) = ^ Üffx) + ^ Üb ( fx).

Therefore any general problem can be solved by applying Eqs. (28) and (30) for all media with the appropriate constraints/boundary conditions.

4.2. Numerical case study and validation

The acoustic transmission analysis framework is applied on the simple macro-scale case study shown in Fig. 4. The system consists of a metamaterial medium made of n =10 unit cells starting at x = 0 till x = nf, where € is the size of the unit cell. The 1D macro-scale coordinate axis and the wave propagation direction are taken along e1 as defined in the unit cell Fig. 1. An acoustic actuation source is applied at x = 0, which prescribes a given displacement at the interface. The

400 800 1200 1600 frequency [Hz]

§ 10^

llO"8 10-12

Horiz. SA Vert.

X DNS X

400 800 1200 1600 frequency [Hz]

Fig. 5. The transmission ratio for applied horizontal (vertical) excitation and measured horizontal (vertical) displacement in the example macro-scale problem using (a) UC2 and (b) UC3 as the LRAM medium computed using the semi-analytical approach (SA) and direct numerical simulation (DNS).

metamaterial medium is bounded on the right side by a semi-infinite homogeneous medium and on the left by an actuation source that applies a prescribed displacement. The impedances of the bounding medium and the actuation source are matched with the impedance of the matrix material within the LRAM. The solution to this problem using the semi-analytical approach is derived in A.1. To verify the semi-analytical solution, a direct numerical simulation (DNS) is carried out using a standard FE software package (COMSOL). The full waveguide structure used in the DNS is built by serially repeating the FE unit cell model used for the dispersion analysis in Section 3.4. Periodic boundary conditions are applied to the top and bottom edges to mimic an infinitely large structure in the vertical direction. An acoustic impedance boundary condition is added at both the ends to account for the impedance of the bounding media and the actuation source. The assembled system has a number of degrees of freedom in the order of 105 .

The transmission analysis using the semi-analytical approach and DNS are first carried out for the LRAM medium consisting of UC2 unit cells. The results are shown in Fig. 5a. The analysis is performed for applied unit horizontal and vertical displacements separately (while fixing the displacement in the other direction). The absolute value of the displacement is measured at the right interface which gives the corresponding transmission coefficient with respect to the applied excitation. The expression for the transmission coefficient obtained using the semi-analytical approach is given by Eq. (A.8). Due to the plane symmetry of the UC2 along el, no hybrid wave solutions exist and a horizontal applied displacement will excite a pure compressive wave and a vertical displacement will excite a pure shear wave. Hence, only the component of the transmission coefficient corresponding to that of the nonzero applied displacement is shown in the figure. The applied frequency range is 350-2100 Hz in order to capture the first two compressive wave bandgaps and the first three shear wave bandgaps (see Fig. 2b). An excellent match between the semi-analytical model and DNS is obtained.

One of the interesting observations is that the shear transmission bandgaps are more pronounced (deeper) than those of the compressive wave. This can be attributed to the fact that the effective shear stiffness of the unit cell is much lower than its compressive stiffness. Due to the requirements of the scale separation, the semi-analytical model is only valid for a relatively stiff matrix material. Beyond this limit, Bragg scattering effects start to play a role, leading to hybridization of Bragg and local resonance effects [37]. This effect does enhance the attenuation performance of the metamaterial but at a cost of overall structural stiffness, hence there is a tradeoff. LRAMs can therefore be employed in applications where a higher structural stiffness are required.

Next, the transmission analysis is performed exactly with UC3 as the LRAM medium (Fig. 5b). Again, the results obtained from the semi-analytical approach match very well with DNS. The UC3 response is similar to UC2 except for the reduced attenuation rate and some additional

°300 400 500 600 700 frequency [Hz]

Fig. 6. The transmission ratio for applied horizontal excitation and measured vertical displacement in the example macro-scale problem using UC3 as the LRAM medium computed using the semi-analytical approach (SA) and direct numerical simulation (DNS).

cross coupling of transmission spectra at the local resonance frequencies. This is due to the fact that a propagating hybrid wave mode will always be excited in UC3 for pure horizontal or vertical excitations at the local resonance frequencies, leading to some residual transmission. The effects due to the hybrid wave modes occurring in UC3 metamaterial are further illustrated in Fig. 6, where the absolute vertical displacement at the right interface is shown, for the unit horizontal applied displacement. The frequency now ranges from 300 to 700 Hz in order to capture the effects due to the first local resonance at 508 Hz. Comparing the results with DNS shows, once again, a perfect match between the two approaches. The peak vertical displacements in the LRAM medium is observed exactly at the first local resonance frequency and drops rapidly away from it, indicating the formation of hybrid modes due to the dynamic anisotropy of the mass density tensor at the resonance frequency. This characteristic response of the UC3 LRAM can lead to interesting applications. As a consequence of horizontal applied displacements on the LRAM, a shear wave will be excited in the homogeneous medium bounding the LRAM. This enables the design of selective mode converters that convert an incident wave mode into another mode upon transmission at particular frequencies for normal wave incidence.

5. Conclusion

A multiscale semi-analytical technique was presented for the acoustic plane wave analysis of (negative) dynamic mass density type LRAMs. It enables an efficient and accurate computation of the dispersion characteristics and acoustic performance of LRAM structures with complex RVE geometries in the low frequency regime. The technique uses a two

scale solution ansatz in which the micro-scale problem is discretized using FE (capturing the detailed micro-dynamic effects of a particular micro-structure) and where an analytical plane wave solution is recovered at the macro-scale analysis. This approach offers two main advantages over traditional Bloch based approaches for dispersion and transmission analysis of LRAMs.

• First, as mentioned earlier, it is computationally remarkably cheap, even in comparison to the reduced order Bloch analysis approaches. The computation of the homogenized enriched material coefficients forms the biggest part of the overall numerical cost which involves meshing, matrix assembly and model reduction. For a particular RVE design, this computation has to be performed only once and therefore constitutes an offline cost. The dispersion spectrum is then computed very cheaply due to the highly reduced nature of the corresponding eigenvalue problem. The mode selection criterion ensures that the resulting model is as compact as possible without sacrificing the accuracy of the solution.

• Second, it permits an intuitive and insightful characterization of LRAMs and its dispersive properties in terms of a compact set of enriched effective material parameters. The only limitation of the method is that it is only applicable towards the analysis of LRAMs in the low frequency regime.

Appendix A

A1. Solution to the acoustic transmission case study problem

The acoustic transmission analysis framework discussed in Section 4 is applied to obtain the results described in Section 4.2 and Fig. 4. The actuation source on the left is indexed as 1, the LRAM medium in the middle as 2 and the homogeneous medium on the right 3. Note that medium 3 has only forward traveling waves since it is semi-infinite. The coordinates of the left (x = 0) and the right interface (x = nf) are represented by 1x and 2x. respectively. Let the prescribed displacement acting on the interface at 1x be denoted as 1 wapp and the

resulting reactive traction be denoted as 1 iapp . The transmitted wave displacement into medium 3, u(2x) as a function of 1 uapp gives the main result of this section.

Applying the boundary conditions on the plane wave solution given by Eq. (30) yields,

. 2 , A -i , 2 ) ,1 .

1 7.1 ,-i _,• „27/2 ),j 1 v\ _2

t,,„ = — im 1 Z • 1 uapp — im 2Z •

(2 f - x) — 2Ib ( - x) - ,

(A.1a)

(A.1b)

The primary result of the work is the derivation of the dispersion eigenvalue problem of the enriched continuum, which gives the fundamental plane wave solution in an infinite LRAM medium. The problem can be conveniently framed in the k - m or m - k form, the latter allowing evanescent wave solutions to be computed. Furthermore, the simple structure of the dispersion model leads to clarifying qualitative and quantitative insights. For instance, the relation between the geometric symmetries, the dynamic anisotropy of the resulting dynamic mass tensor and the nature of the dispersive behavior becomes clear, i.e. the formation of selective/total bandgaps and/or hybrid wave modes at certain frequencies. Several case studies were used to illustrate the approach. The effective continuum description was combined with a general plane wave expansion resulting in a modified transfer matrix method on enriched media. The analysis was restricted to normal wave incidence, but it can easily be extended to the general case of an oblique incidence.

As discussed earlier, complex designs induce extended micro-dynamics, which manifests itself at the macro-scale in two important ways. First, in the proliferation of local resonance modes, which trigger the formation of additional bandgaps in the dispersion spectrum. And second, in the dynamic anisotropy of the material response which leads to the formation of selective bandgaps and hybrid wave modes. These effects can be exploited towards developing more advanced filters and frequency multiplexers, which target specific wave modes at multiple desired frequencies. The developed framework was applied to case studies revealing interesting phenomena such as pronounced shear transmission bandgaps compared to the compressive ones (due to the difference in the elastic stiffness). The localized hybrid mode formation in considered LRAM unit cell was used to demonstrate a potential application towards selective mode conversion, where the energy of the incident normal wave mode can be channeled into a different mode within a certain frequency range. Hence, the presented semi-analytical approach serves as a valuable tool for optimization and rapid prototyping of LRAMs for specific engineering applications.

Acknowledgments

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no [339392].

I( 2 x) = 3Mf (2 x) = 2 Mf (2 x) + 2 Mb (2 x) ,

0 = im2Z • (2wf(2x) — 2Ib(2x)) — im 3Z • fx).

(A.1c)

(A.1d)

Note that the overbar is not applied on medium 1 and 3 to indicate that they are not homogenized quantities. Solving for 2 «b (2 x) with respect to 2 f2 x) using Eqs. (A.1c) and (A.1d) gives

2 Ib (- x) = ARe • 2 f2 x) ,

where,

Are = (2 Z + 3 Z)-1 • (2 Z — 3 Z)

is the effective reflection coefficient tensor at the given interface. Applying the transfer operation as defined by Eq. (28) gives 2 f2 x) = 2T(«0 • 2 Mf( 1 x) and 2 Mb( 1 x) = 2T(n€) ■ 2 ub(2 x). Combining the result with Eq. (A.2) gives

2 Ib( - x) = 2 T(nf) • ARe • 2 T(«0 • 2 f1 x) .

Substituting the above expression into Eq. (A.1a) and solving for 2 Uf(1 x)

2 Uf(1 x) = (i + 2T(«^) • ÂRe • 2 W• 'Û^app • (A.5)

Now, applying the transfer operation to Eq. (A.5), the solution 2 3f(2 x) is obtained as

2 f2 x) = 2 T(«€) • (i + 2 T(«0 • ARe •2 T(«0 - 1 • 1 u

Finally, using Eqs. (A.2) and (A.1c) and the fact that f3 x) = w(2 x) gives the desired result

îf (- x) = (i + ARe - • 2 T(nf) • (i + 2 T(nf) • ARe • 2 T( nf) - • 1

: ATr • M app

where ATr is the effective transmission coefficient between the applied displacement and the transmitted wave in medium 3, expressed as

ATr = (l + ARe ) • 2 T( n€) • (1 + 2 T(nO • ARe • 2 T(n^) .

References

[1] Deymier PA. Acoustic metamaterials and phononic crystals, 12. Springer Series in solid-state sciences; 2013 .

[2] Haberman MR, Guild MD . Acoustic metamaterials. Phys Today 2016;69:42-8.

[3] Cummer SA, Christensen J, Alu A. Controlling sound with acoustic metamaterials. Nat Rev Mater 2016;1(16001):528-51 .

[4] Ma G, Sheng P . Acoustic metamaterials: from local resonances to broad horizons. Sci Adv 2016;2(2) e1501595.

[5] Liu Z, Zhang X, Mao Y, Zhu YY, Yang Z, Chan CT, et al. Locally resonant sonic materials. Science 2000;289:1734-6,

[6] Fang N , Xi D, Xu J, Ambati M, Srituravanich W, Sun C, et al. Ultrasonic metamaterials with negative modulus. Nat Mater 2006;5:452-6.

[7] Larabi H, Pennec Y, Djafari-Rouhani B, Vasseur JO , Multicoaxial cylindrical inclusions in locally resonant phononic crystals. Phys Rev E 2007;75(6) ,

[8] Gazalet J, Dupont S, Kastelik JC, Rolland Q, Djafari-Rouhani B , A tutorial survey on waves propagating in periodic media: electronic, photonic and phononic crystals. perception of the bloch theorem in both real and fourier domains. Wave Motion 2013;50:619-54.

[9] Kafesaki M, Economou EN. Multiple-scattering theory for three-dimensional periodic acoustic composites. Phys Rev B 1999;60(17):1993-12001.

[10] Liu Z, Chan CT, Sheng P. Elastic wave scattering by periodic structures of spherical objects: theory and experiment. Phys Rev B 2000;62(4):2446-57.

[11] Mead DJ. A general theory of harmonic wave propagation in linear periodic systems with multiple coupling. J Sound Vib 1973;27(2):235-60.

[12] Collet M, Ouisse M, Ruzzene M, Ichchou MN , Floquet-bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems. Int J Solids Struct 2011;48(20):2837-48.

[13] Kulpe JA, Sabra KG, Leamy MJ. Bloch-wave expansion technique for predicting wave reflection and transmission in two-dimensional phononic crystals. J Acoust Soc Am 2014;135(4):1808-19.

[14] Farzbod F, Leamy MJ , Analysis of bloch ' s method and the propagation technique in periodic structures. J Vib Acoust 2011;133(3) 031010.

[15] Zhou X, Liu X, Hu G. Elastic metamaterials with local resonances: an overview. Theor Appl Mech Lett 2012;2(4) 041001.

[16] Sheng P. Introduction to wave scattering, localization, and mesoscopic phenomena. Academic Press; 1995.

[17] Wang J, Michelitsch TM , Gao H, Levin VM. On the solution of the dynamic eshelby problem for inclusions of various shapes. Int J Solids Struct 2005;42:353-63.

[18] Auriault J-L, Boutin C , Long wavelength inner-resonance cut-off frequencies in elastic composite materials. Int J Solids Struct 2012;49(23-24):3269-81.

[19] Boutin C Hans S3 Chesnais С Generalized beams and continua. dynamics of reticulated structures. Mech Gen Contin 2010,;21:131-41

[20] Pham K Kouznetsova VG3 Geers MGD . Transient computational homogenization for heterogeneous materials under dynamic excitation. J Mech Phys Solids 2013;61(11):2125-46■

[21] Sridhar A, Kouznetsova VG3 Geers MGD • Homogenization of locally resonant acoustic metamaterials: towards an emergent enriched continuum. Comput Mech 2016;57(3):423-35 .

[22] Sheng P3 Mei J3 Liu Z, Wen W • Dynamic mass density and acoustic metamaterials. Phys B: Condens Matter 2007;394(2):256-61

[23] Lai Y3 Wu ^ Sheng P3 Zhang Z • Hybrid elastic solids. Nat Mater 2011;10:620-4•

[24] Michelitsch TM3 Gao H, Levin VM • Dynamic eshelby tensor and potentials for ellipsoidal inclusions. Proc R Soc, A 2003;459:863-90•

[25] Eshelby JD . The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc R Soc A 1957;241(1226):376-96.

[26] Reda H Ganghoffer JF3 Lakiss H • Micropolar dissipative models for the analysis of 2d dispersive waves in periodic lattices. J Sound Vib 2017;392:325-45.

[27] Reda H, Rahali Y, Ganghoffer JF, Lakiss H. Nonlinear dynamical analysis of 3d textiles based on second order gradient homogenized media. Comp Struct 2016;154:538-55 .

[28] Eringen AC. Microcontinuum field theories. Springer-Verlag New York; 1999.

[29] Kouznetsova VG3 Brekelmans WAM Baaijens FPT . An approach to micro-macro modeling of heterogeneous materials. Comput Mech 2001;27(1):37-48.

[30] Nemat-Nasser S, Hori M . Micromechanics: overall properties of heterogeneous materials. North-Holland series in applied mathematics and mechanics Elsevier; 1993.

[31] Milton GW3 Willis JR • On modifications of newton 3 second l aw and linear continuum elastodynamics. Proc R Soc 2007;463(2079):855-80.

[32] Bückmann T Kadic M, Schittny R3 Wegener M • Mechanical metamaterials with anisotropic and negative effective mass-density tensor made from one constituent material. Phys Status Solidi (B) 2015;252(7):1671-4.

[33] Nassar H He Q3 Auffray N • Willis elastodynamic homogenization theory revisited for periodic media. J Mech Phys Solids 2015;77:158-78.

[34] Hussein MI . Reduced bloch mode expansion for periodic media band structure calculations. Proc R Soc A 2009;465(2109):2825-48.

[35] Krattiger D, Hussein MI. Bloch mode synthesis: ultrafast methodology for elastic band-structure calculations. Phys Rev E 2014;90(6) 063306.

[36] Lowe MJS . Matrix techniques for modeling ultrasonic waves in multilayered media. IEEE Trans Ultrason, Ferroelectr Frequency Control 1995;42(4):525-42•

[37] Yuan B , Humphrey VF3 Wen J3 Wen X • On the coupling of resonance and bragg scattering effects in three-dimensional locally resonant sonic materials. Ultrasonics 2013;53(7):1332-43.