Homogenization of locally resonant acoustic metamaterials towards an emergent enriched continuum.

This contribution presents a novel homogenization technique for modeling heterogeneous materials with micro-inertia effects such as locally resonant acoustic metamaterials. Linear elastodynamics is used to model the micro and macro scale problems and an extended first order Computational Homogenization framework is used to establish the coupling. Craig Bampton Mode Synthesis is then applied to solve and eliminate the microscale problem, resulting in a compact closed form description of the microdynamics that accurately captures the Local Resonance phenomena. The resulting equations represent an enriched continuum in which additional kinematic degrees of freedom emerge to account for Local Resonance effects which would otherwise be absent in a classical continuum. Such an approach retains the accuracy and robustness offered by a standard Computational Homogenization implementation, whereby the problem and the computational time are reduced to the on-line solution of one scale only.

Introduction

The study of the elastodynamics of heterogeneous materials has led to the discovery of a special class of a composite material known as Acoustic Metamaterials [1]. They exhibit remarkable acoustic properties ranging from near zero transmissibility [2], enhanced absorption [3], negative dynamic mass density [4]/ bulk modulus [5], negative refractive index [6], super anisotropy, zero rigidity [7] etc. These exotic phenomena have numerous potential applications such as low frequency noise attenuation [3], isolation of civil structures from seismic waves [8], superlenses with a resolution beyond the Rayleigh limit [6, 9], waveguides that can be used to channel acoustic waves, etc.

Two important physical phenomena are responsible for the extraordinary properties of Acoustic Metamaterials, Local Resonance and Bragg Scattering [10]. These two phenomena operate at different length scales, whereby Bragg scattering is dominant for wavelengths (of the propagating wave) of the same order as the size of microstructural phases and Local Resonance at larger wavelengths. This work exclusively deals with the modeling of the latter phenomena, which are typically applied in the lower frequency regime. The subclass of Acoustic Metamaterials exhibiting Local Resonance is known as locally resonant acoustic metamaterials (LRAM). A typical unit cell microstructure of a LRAM consists of a matrix with an embedded inclusion or a substructure [10]. These inclusions/substructures consist of two parts, a central region with high mass density supported by a surrounding highly compliant region (for eg. rubber coated lead inclusions [2], also see Fig. 4). This enables the unit cell to exhibit low frequency localized vibration modes (see Fig. 5) that strongly couple to the long wavelength propagating wave in the matrix at the resonance frequency. The strong coupling around this frequency is what is responsible for the Local Resonance phenomena [2, 4].

Fig. 4

Example of a LRAM unit cell designed by [2]

Fig. 5

Local resonance mode shapes and their associated eigenfrequencies and momentum coupling coefficients of the unit cell under consideration

A plethora of techniques is available for modeling the elastodynamics of heterogeneous materials although not all are equally suitable for describing LRAM. Direct numerical simulation (DNS) of LRAM structures using finite element method (FEM) is highly unpractical due to the large scale difference involved in such problems. Therefore it is necessary to develop efficient methods with more ingenuity. The Bloch-Floquet (or simply Bloch) theory [10] provides a general solution for the propagation of free waves (i.e. steady state waves in an infinite medium) in any periodic heterogeneous material. Substitution of the Bloch solution into the governing equations reduces the entire analysis to a single unit cell. The unit cell problem can then be solved numerically using various discretization methods such as Plane Wave Expansion [11], Variational Method [12], Finite Difference Time Domain [1] and FEM [13]. FEM provides a high flexibility in the design of the unit cell topologies at the cost of convergence with respect to the mesh size. This limitation can be overcome to a great extent by using model reduction techniques [14]. A homogenization scheme based on the Bloch theory was proposed by Willis [15, 16] and revisited again in [17]. The techniques based on Bloch theory has been highly successful in the study of free wave propagation but is mostly limited to this case. It is difficult if not impossible to account for complex transient loading and macroscopic boundary effects using this theory. Apart from Bloch theory, another generalized solution for transmission of free waves in a heterogeneous medium is given by the Multiple Scattering Theory [18]. It can be used to predict macroscopic transmission spectra and provides superior convergence with respect to discretization size, yet it is highly restricted to simple unit cell topologies. To the best of our knowledge it has only been implemented for spherical inclusions and granular media.

Another modeling approach is provided by Enriched or Micromorphic Continuum Theory [19], first proposed for elastodynamics by Mindlin [20]. It introduces additional macroscopic kinematic fields that account for the internal microscale dynamics in an otherwise homogeneous macroscopic medium. A new material property called ‘micro-inertia’ is postulated that characterizes the microdynamics. However, the only notable attempt at developing an enriched model that is especially capable of accounting for Local Resonance has been made by Sun et al. [21]. Homogenization, namely Asymptotic Homogenization [22-28] has also been a highly successful approach for modeling dispersion behavior in heterogeneous materials but few works exist that are specifically suitable for modeling LRAM.

The techniques developed thus far have been successful in studying LRAM behavior but do not provide a generalized description. A more comprehensive modeling technique for LRAM should enable the description and analysis of complex microstructure topologies, transient response at both scales and finite macroscopic structures with various boundary conditions. Motivated by this challenge, this contribution presents a novel framework for modeling LRAM that is not only general in its description but also highly efficient, enabling a fast and straightforward numerical implementation.

An extended first order Computational Homogenization framework [29] is taken as the point of departure to setup the multiscale transient dynamic problem. Except for linear elasticity and relaxed scale separation assumption the framework is general and the full balance of the linear momentum is solved at both scales.The relaxed scale separation principle introduced here retains the long wavelength (quasistatic) assumption on the matrix material but relaxes it on the inclusions. This accounts for the transient dynamic behavior of the microstructure characterizing micro-inertia effects, especially Local Resonance. A Computational Homogenization framework can be efficiently combined with FEM techniques for multiscale problems, which in turn gives the freedom to incorporate complex microsctructure topologies and finite macrostructure geometries, arbitrary transient excitation and sophisticated boundary conditions. The first use of Computational Homogenization to model LRAMs was made by Pham et al. [30]. The approach proposed in this work is distinctly different as it aims at obtaining a closed-form description of the macroscopic continuum, which is enriched to incorporate the effect of microscale dynamics. To eliminate the on-line (expensive) solutions of the microscale problems at each time step, typically used in a computational homogenization approach, a technique called the Craig Bampton Mode Synthesis [31] is employed. It entails a decomposition of the solution into two parts, the quasistatic response and internal dynamics, which makes it possible to condense the behavior of the microscale model upto the macroscopic level by applying homogenization. The method employs eigenmodes to extract the relevant dynamics of the microscale problem, which effectively captures Local Resonance effects with a minimum set of degrees of freedom. A highly compact closed-form description of a linear elastic LRAM results and the corresponding expressions effectively represent an enriched continuum where the emerging additional field accounts for the microscopic Local Resonance phenomena. This approach therefore provides a simple, efficient and general description of a linear elastic LRAM. This not only enables an efficient numerical implementation but also provides an intuitive understanding of the behavior of such materials.

The formulation used here to setup the micro and macroscale problems is defined for a Cauchy continuum in a two-dimensional space. The method can be easily extended to three-dimensions or other specialized continua such as shells, beams etc. The paper is organized as follows. Section 2 presents the overall methodology, including the homogenization framework for the problem and the techniques used to reduce the model towards an enriched continuum. In Sect. 3, the model is numerically validated against DNS for a 1D compressional wave test on a well known example of a LRAM structure. The conclusions are given in Sect. 4.

The following notations are used throughout the paper to represent different quantities and operations. The Cartesian basis vectors are given by . Unless otherwise stated scalars, vectors, second- and fourth-order Cartesian tensors are generally denoted by a (or A), , and respectively. The standard tensor operations are denoted as follows, dyadic product: , dot product: and double contraction: (Einstein summation is used here and for all tensor operations). The conjugate of any second order tensor, is indicated as . Matrices of any type of quantity are in general denoted by and the special case of a column matrix is denoted by . Submatrices of matrix or are denoted by the left superscript and respectively. Transpose of a matrix is given as .

Homogenization and reduction methodology: towards an enriched continuum

The classical first order homogenization framework extended to the transient dynamic case is used to setup the multiscale problem. The full balance of linear momentum is considered at both scales. Scale transition relations are formulated that dictate the coupling between the scales. Once the framework is defined, the focus is shifted to the microscale problem which is reformulated in a discretized form. Using the Craig Bampton technique, a compact reduced model of the microstructure is obtained. The scale transition relations are then applied to upscale the reduced microscale balance equations, yielding the governing macroscale continuum equations, providing a closed form description of the enriched macroscopic continuum.

Homogenization framework

Separation of scales

In the classical first order homogenization scheme, the separation of scales principle requires the following relation to hold true for a heterogeneous material with n microstructural constituents (cavities, grains, inclusions, matrix etc.):where is the typical size of the jth microstructure constituent and the corresponding shortest characteristic wavelength of the microstructural constituent for a given applied excitation. Under this assumption, the micro inertial response of the microstructure becomes negligible leading to a purely quasistatic response. Therefore a more relaxed scale separation principle is adopted here. Let and be the number of microstructural phases constituting the heterogeneities and the core matrix, respectively:where and are the shortest characteristic wavelengths in the jth and kth constituents of the matrix and the heterogeneity for a given applied excitation, respectively. The long wavelength approximation still applies to the matrix whereas a more relaxed hypothesis holds for the heterogeneities. The microstructural lengths can now scale with the wavelengths associated to heterogeneities, incorporating possible micro-inertia effects.

Macroscale problem

Let represent the domain of the macroscale problem and its boundary. Let give the position vector of any point in this domain. A Cauchy continuum is assumed with the displacement vector and its gradient representing the primary kinematic fields. It will be shown later that static rotational equilibrium is not necessarily a priori satisfied at each continuum point which requires the use of the full displacement gradient instead of the symmetric linear strain tensor. The governing equation in the absence of external body forces is given by the linear balance of momentum,where and are the macroscopic stress tensor and momentum respectively. Appropriate initial and boundary conditions should be applied but these are not explicitly stated here since they are not important for the subsequent derivations. The macroscale constitutive response1 is derived from the solution of the microscopic problem via homogenization. In classical Computational Homogenization schemes, this macroscopic constitutive response is not obtained through a closed-form equation. This is different (and better) for the present case, see Sect. 2.3, which constitutes an important step forward.

Microscale problem

To each material point of , a fine scale domain with boundary is associated, where the microscale problem is defined. The domain is selected such that it captures the local microstructural effects at that point. It is termed the representative volume element (RVE). For (locally) periodic microstructures, the RVE is defined as the unit cell that spans the (local) microstructure. For the sake of simplicity, no specific subscripts are used to indicate the variables associated to the microstructural domain. Let the total volume and the infinitesimal volume element of be represented by V and respectively and an infinitesimal surface element of be represented by . A Cauchy continuum is assumed for the microscopic problem with the kinematic field variables given by displacement and linear strain . In order to capture the micro inertial effects, the full balance of momentum has to be considered in the RVE,A linear elastic material is considered. For every material constituent domain , the following constitutive relations hold, Where and and stand for the mass density and the elastic material stiffness of the constituent respectively. A perfect bonding condition is assumed; Here are the unit normal vectors to the interface between both constituents. The initial and boundary conditions on the RVE required to solve the micro problem will be defined upon introducing the scale transition relations.

General microscopic kinematics

The kinematics of the RVE corresponding to a point of is given by the following first order representation of the microscopic kinematics at that point,Here, is a reference vector whose definition will be given later and , called the microfluctuation field, represents the fine scale variations due to the microstructure heterogeneities. This field provides the necessary kinematic degrees of freedom to describe the Local Resonance phenomena.

Scale transition relations

With the introduction of the microfluctuation field , additional constraints in form of boundary conditions on the RVE are required to ensure the well posedness of the problem. These conditions usually follow from the (chosen) relations coupling both scales. They essentially represent the kinematic coupling between the macro and microscale and are called downscaling relations. In the micro-macro direction upscaling relations recover the macroscopic stress and momentum from the solution of the RVE boundary value problem. They result by inserting the downscaling relations into the Hill-Mandel macrohomogeneity condition, generalized to the transient dynamic case. The up and downscaling relations together constitute the scale transition relations that dictate the coupling between the two scales.

Downscaling relations (kinematic boundary conditions): The condition on the macroscopic displacement is formulated as the overall rigid body displacement of the RVE. This is equivalent to constraining the microfluctuation at a single arbitrary point, say , on the RVE boundary to zero.The second condition on the macroscopic displacement gradient follows the established averaging theorem [32]Substituting the RVE kinematics given by Eq. (7) into the above expression and simplifying results inThis gives the constraints on the microfluctuation field, which is easily converted to a boundary integral by applying Gauss theoremwhere is the outward normal to the RVE boundary. Relation (11) gives the minimal kinematic constraint that enforces Eq. (9) for a given macroscopic displacement gradient. It has been shown in the literature that this constraint is too weak and generally leads to poor results for small RVEs. Therefore a stronger form known as periodic boundary condition (pbc) is applied instead, which provides much better RVE size convergence [33]. The assumption of periodic boundary displacements can be justified for a transient problem if the matrix is assumed to behave quasi-statically. Indeed, this holds true in the case of the relaxed scale separation principle given by Eq. (2), which retains the long wavelength assumption on the matrix. The periodic boundary conditions can be formulated as follows. Considering a RVE with a periodic rectangular shape (in 2D), the boundary is split into four edges: left, right, bottom and top, denoted by L, R, B and T, respectively. The four vertices are denoted by , , and (see Fig. 1). Let the position vectors of the vertices be represented by , , and , respectively. The normal unit vector at every corresponding pair of points on opposite boundary sides satisfies, and . By constraining the microfluctuations on the edges to be periodic, i.e.

Fig. 1

Sketch of the boundaries of a RVE and its normal vectors

Equation (11) is automatically satisfied. Making the choice, in Eq. (8) and substituting it in Eq. (12), gives the following condition on the microfluctuations at nodes p, p and p (which will be used in a later section)Upscaling relations (Hill-Mandel condition): According to the principle of virtual work, the total internal virtual work must be conserved for a dynamic system for any imposed kinematics. The Hill-Mandel principle extends this concept to relate the two scales where the volume averaged virtual work of the RVE is equated to its macroscopic equivalent in a material point. Therefore, equating the volume average of the virtual work of Eqs. (3) and (4) givesRewriting the left hand side of the above expression in terms of the virtual work of the external boundary tractions on the RVE gives,Substituting Eq. (7) into the above expression and making use of Eqs. (12) and (13) givesWhich upon simplification reads,This gives the corresponding upscaling relations. Note that in homogenization of static problems, vanishes from Eq. (17) due to static equilibrium (zero net traction on the boundary) of the RVE. This is no longer the case in a dynamic setting where the net tractions do not vanish, see Eq. (18), indicating the dependency of the macroscopic stress on and also on the micro-inertial effects. By employing Gauss theorem and the microscopic balance given by Eq. (4), the boundary integral form of Eqs. (17) and (18) can be written in the equivalent volume integral form as follows,Note that Eq. (19) highlights the explicit coupling of the macroscopic stress to the microscopic momentum, again indicating the influence of micro-inertial effects on the macroscopic stress.

Discretization and model reduction

Using standard discretization techniques, e.g. based on the finite element method (FEM), the discrete balance of momentum of the RVE can be written as followsHere and represent the nodal displacements, accelerations and applied forces respectively in which the quantities associated to each node are represented by a vector and assembled in a vector column fashion. represents the stiffness tensor matrix and the mass tensor2 matrix which transform the nodal displacement and acceleration vectors respectively into internal nodal force vectors. In the derivations that ensue, standard Galerkin projection is used to perform model reduction.

Periodic boundary condition

In order to apply the periodic boundary conditions it is necessary to partition the system into tied (constrained) and retained nodes denoted by the characters ‘t’ and ‘r’ respectively. The tied nodes constitute the right and the top edge nodes along with vector point and the retained nodes constitute the remaining boundary nodes and the nodes in the RVE interior. The discrete form of Eq. (12) can be written in terms of displacements as follows,where denotes a column with unit scalar entries of the same size as and . The above expressions give the constraint relations on the tied nodes in terms of the retained nodes. Let represent the scalar reduction matrix reflecting this linear transformation;Applying the reduction (25) on Eq. (21) leads to the reduced discrete governing equations in terms of the retained nodeswhere

Craig Bampton reduction

The Craig Bampton Mode Synthesis [31] (or Craig Bampton for/CBMS short) is a popular substructuring technique used in structural dynamics to obtain reduced models of complex assembled systems (e.g. cars, planes etc). It involves a reduced description of the dynamic response of the interior of every subsystem or substructure with respect to prescribed displacements at the external boundary of each subsystem. The total response is then obtained by assembling the individual reduced substructure models at the boundaries separating the substructures. This concept can be applied to a multiscale framework where the microscale problem is treated as a substructure of the macroscale problem. Eigenmodes of the interior dynamics are used as a reduced basis, which perfectly capture the Local Resonance phenomena.

With the periodic boundary conditions incorporated, the Craig Bampton procedure can now be applied. The total dynamic response of any structure under prescribed displacements (velocities, accelerations) can be expressed as a superposition of its corresponding quasistatic response and its internal dynamics spanned by a minimal set of eigenmodes of the structure with the prescribed nodes fixed. This is illustrated in Fig. 2. This decomposition is essential since it separates the Local Resonance effects from the rest of the mechanical response thereby enabling the construction of a reduced model. The quasistatic response gives the instantaneous mechanical response whereas the internal dynamics represents the inertial response due to Local Resonance. The two subproblems are next treated independently and later superposed to obtain the total solution. Note that the incorporation of internal dynamics is made possible due to the assumption of the relaxed scale separation principle given by Eq. (2). In general, when applied incorrectly, the Craig Bampton reduction could lead to stiff and sometimes incorrect responses since the eigenmodes are computed on a constrained system, which may not capture the actual dynamics. However this is not an issue for homogenization as long as the prescribed displacement nodes lie on the matrix (which does not exhibit relevant dynamical behavior in the frequency ranges satisfying the relaxed separation of scales) and not the inclusion.

Fig. 2

Illustration of the Craig Bampton decomposition. The total dynamic response of an RVE to prescribed boundary displacements can be represented as a superposition of its quasistatic response and its internal dynamics spanned by a set of eigenmodes with the prescribed nodes fixed

In order to proceed with the derivations, the retained nodes must be further partitioned into prescribed and free nodes denoted by ‘p’ and ‘f’, respectively. The prescribed nodes include nodes p, p and p and the free nodes consist of the remaining nodes. The partitioned form of Eq. (26) is written aswhere represents a column of length (where indicates the number of free nodes) in which each entry is a zero vector. It can be seen that no external forces act on the free nodes.

Quasistatic Response: The quasistatic response is recovered by omitting the mass contribution in Eq. (28), enabling to solve for the free degrees of freedom using a static condensation procedure. The solution can then be represented in terms of the prescribed displacements as follows (subscript ‘qs’ is used for quasistatic)The term represents the static condensate (also known as the Schur complement) of the stiffness matrix , and represents a matrix with each entry containing a second order identity tensor. The expression for readsUsing Eq. (29) as a reduced basis, model reduction of Equation (28) gives the governing equations for the quasistatic response expressed in displacements of the prescribed nodes onlywhere,Here, the superscript represents the transpose of a matrix in which each of the tensor entries are conjugated.

Internal dynamics: The internal dynamics is added to the solution in order to compensate for neglecting the contribution of the inertial forces when calculating the quasistatic response. This correction is essential in accounting for Local Resonance. The solution of the dynamic subproblem is expressed as a superposition of eigenmodes (here , being the number of free degrees of freedom), computed with prescribed nodes fixed. The expression for this readswhere represents the vector matrix containing the eigenmodes, the corresponding column of generalized scalar displacements, i.e. the amplitude of the eigenmodes also referred to as modal displacements and represents a zero vector matrix. The eigenvalue problem and the mass normalization condition of each eigenmode , where are respectively, given aswhere is the eigenfrequency corresponding to and represents a zero vector column of length . To yield a computationally efficient model, the basis should be constructed such that it only contains the essential (or excited) modes that trigger Local Resonance. A more quantitative mode selection criteria is briefly described in Sect. 2.3. In addition, the retained eigenmodes must also satisfy the scale separation principle given by Eq. (2). Applying Eq. (33) to reduce Eq. (21) with the help of Eqs. (34) and (35) giveswhere is a zero scalar column of length and,Note that Eq. (36) describes a set of uncoupled spring mass systems with eigenfrequencies given by . This gives the simplest description of Local Resonance without any coupling to the macroscopic dynamics. This coupling will be obtained by combining the solutions of the two subproblems.

Linear superposition: The full solution is now expressed as a superposition of the quasistatic response and the internal dynamics . Combining Eqs. (29) and (33) gives,Projecting Eq. (28) on the subspace of Eq. (38) gives the reduced coupled dynamic model of the RVE.whereis the vector matrix providing the coupling between the two equations.

Emerging enriched continuum

The discrete form of the downscaling relations given by Eq. (13) with account for Eq. (7) is represented aswhere is an identity column of length 3 and , where is a unit scalar column of size 3. Similarly, the discrete form of the upscaling relations Eqs. (17) and (18) is given asThe macroscopic constitutive relations can now be recovered by substituting Eqs. (40) and (42) into the above expressions:Here stands for the zero tensor and stands for the left conjugate of a higher order tensor defined as . As shown in the above expression, several terms are zero or can be neglected. The corresponding justification for each of these terms is given below:

These terms are identically zero because is in fact the null space of as it describes the rigid body modes imposed by the prescribed nodes.

This term gives the elastic inertia (inertial contribution to deformation modes) and for the materials of concern here it is negligible compared to the elastic stiffness (i.e. ). Although the total term, can become significant at higher applied frequencies, (since ), these frequencies approach the homogenizability limit (where scale separation is violated), hence it is justified to drop this term in the considered regime.

These are the cross coupling terms between and and and . They are in general not zero and depend on the choice of , which was yet unspecified. The choice of does not influence the accuracy of the final solution. Hence, the value of is simply chosen such that these terms become zero. The physical interpretation of these terms will be unraveled in future work.

The homogenized constitutive expressions given by Eqs. (45) and (46) can now be written in a compact closed-form as follows,where, Here represents the effective elasticity tensor, is the average mass density of the RVE. The terms and represent the coupling of the sth modal degree of freedom to the macroscopic stress and momentum respectively. They can be used to determine the relevant eigenmodes for the reduced basis. Local resonance modes must have at least one nonzero (and finite) coefficient in or . The displacement gradient in Eq. (47) can be replaced with the linear strain tensor due to the symmetries present in the homogenized elasticity tensor . However can be asymmetric since rotational eigenmodes can exist, introducing angular inertia into the system violating the rotational equilibrium condition accompanying the symmetry of the stress tensor. Hence, this triggers the asymmetry of the macroscopic stress tensor in materials with micro-inertia, in general. In the example problem considered in this work, the rotational modes have a negligible Local Resonance amplitude and thus the effect of angular inertia does not contribute here. They will be discussed in more detail in forthcoming work.

Finally, the balance of momentum of the internal dynamics needs to be expressed in its homogenized continuum notations. Applying Eqs. (42) to (39) givesThe emerging macroscopic governing equations of the enriched continuum can be summarized below

Macroscopic balance of momentum. Microscopic balance of momentum. (To be solved at the macroscale)Homogenized constitutive relations.

Note that is an emergent field variable, which effectively ‘enriches’ the macroscopic continuum with the micro-inertia effects, in the micromorphic sense as initially defined by Eringen [19]. It constitutes an internal field variable that does not explicitly appear in the classical macroscopic balance equation, and is therefore appropriately termed as internal dynamics.

The entire procedure can be summarized as follows. Starting from the discretized balance of momentum of the RVE, periodic boundary conditions are first applied to the boundary nodes. The Craig Bampton decomposition is introduced by expressing the solution as the superposition of the quasistatic response and the internal dynamics represented by an optimum eigenmode basis. This expression is then used to obtain a reduced balance of momentum of the system. The scale transition relations are then applied to transform the reduced balance of momentum into the homogenized enriched macroscopic continuum equations. The eigenmode problem is typically solved off-line, in advance. The macroscale on-line solution procedure therefore reduces to a single scale enriched problem. The procedure is illustrated in Fig. 3.

Fig. 3

Summary of the main steps involved in computing the enriched continuum model

Numerical validation

Local resonance acoustic metamaterial (LRAM)

For the purpose of validation of the presented model, a well known and understood LRAM design of [2] is adopted. The structure of the unit cell and the macroscopic lattice structure are shown in Fig. 4. It consists of an epoxy matrix with an embedded lead inclusion coated with soft rubber. In [2] spherical inclusions were used. Since derivations are carried out in 2D space, the inclusions are modeled as infinite cylinders instead. The geometric and material properties of the LRAM are given in Table 1. The high compliance of rubber combined with the high mass density of lead compared to the epoxy triggers the observed low frequency localized resonance modes.

Table 1

The geometric and material properties of the unit cell under study

Geometrical properties		Material properties
Radius of lead inclusion	5 mm	Parameter	Lead	Rubber	Epoxy
Thickness of rubber coating	2.5 mm	ρ(×10^3kg/m^3)	11.6	1.3	1.18
Length of the unit cell	20.1 mm	λ (N/m^2)	4.23×10^10	6×10^5	4.43×10^9
		μ (N/m^2)	1.49×10^10	4×10^4	1.59×10^9

The FEM model of the RVE was constructed using the Marc Mentat software package. The discretized system contains close to 6000 degrees of freedom. An eigenvalue analysis was carried out on the system and the first 10 in plane eigenmodes were extracted. Out of these 10 modes, only 2, 3, 5 an 6 are Local Resonance modes with finite values of the momentum coupling coefficient. The remaining modes have either zero or negligible values of their respective stress and momentum coupling coefficients. Modes 2 and 3 and 5 and 6 are degenerate with eigenfrequencies 360 Hz and 1239 Hz respectively. The mode shapes, the eigenfrequencies and the momentum coupling coefficients of each mode are given in Fig. 5. Modes 3 and 5 are coupled along the direction and modes 2 and 6 are coupled along the direction. The mode shapes overlap strongly with those found in [2], hence it can be concluded that these are indeed the Local Resonance modes.

A simple estimation can be made to verify if the RVE under consideration satisfies the scale separation criteria. By neglecting the Local Resonance effects, an approximate estimation of the matrix resonance frequency readswhere L is the length of the RVE and is the speed of sound for a given wave mode ( indicates compressive wave and , shear wave, etc). The above expression gives the order of the frequency at which the matrix starts behaving dynamically and the relaxed scale separation is violated. Therefore, it can be used to obtain an estimation of the limiting frequency. Substituting the values of the respective quantities into the above expression gives a value of the order of Hz for both shear and compressive wave modes. The frequencies of interest determined by the Local Resonance eigenfrequencies (360 and 1239 Hz) lie well within the limiting frequency. Thus the LRAM under consideration is well suited for applying the proposed homogenization method.

Macroscopic problem construction

The developed homogenized enriched continuum is next verified against DNS. For the sake of simplicity, the macroscopic problem is restricted to a 1D compressional wave propagation test. The DNS model is constructed by sequentially stacking 100 RVEs as shown in Fig. 6a. In order to ensure an effective 1D macroscopic behavior for this 2D system, periodic boundary conditions are applied to the top and bottom edges of the RVEs to mimic an infinitely large structure in the vertical direction. The right edge of the system is fully clamped, and prescribed displacements are applied on the left edge. A harmonic loading function of frequency is applied on the left end for 2 time periods as shown in Fig. 6b. The finite element model of the total DNS system contains over 500,000 degrees of freedom.

Fig. 6

Construction of the macroscopic compression wave test using DNS and the enriched homogenized model. The right edge is constrained and a prescribed displacement is applied on the left edge. a Macroscopic problem setup. b Input displacement function

The homogenized model of the system is constructed by discretizing a 1D enriched continuum using 100 linear finite elements with 101 nodes. The problem is discretized in time using an implicit Newmark scheme. The nodal degrees of freedom consist of the horizontal macroscopic displacement and the two generalized modal displacements associated with Mode 3 and 5 (only the modes providing coupling in direction are retained since the problem considered here only involves wave propagation in this direction). The right most node is constrained and the prescribed displacements are applied on the left most node (see Fig. 6a).

Four transient simulations tests were performed at different excitation frequencies , which are selected to probe various regions of interest. It is well known in the literature [2], that the Local Resonance coupling is strongest in the frequency regions above each eigenfrequency known as ‘stop bands’. For the simulation, the first excitation frequency is selected at 200 Hz which is far below the first eigenfrequency where the Local Resonance effects should have a weak influence. The response is expected to be predominantly quasistatic at this frequency. The second frequency is selected at 450 Hz which is in the first stop band. The Local Resonance coupling will be strong there. The third frequency is selected at 800 Hz which is in between the two stop bands and the fourth is selected at 1330 Hz in the second stop band. The results are shown in Fig. 7. Two plots of the horizontal displacement are displayed for each excitation frequency, one as a function of position after two time periods and the second as a function of time at x m. The homogenized solution matches the DNS results in an excellent manner, thus validating the developed method.

Fig. 7

Plot of horizontal nodal displacement [m] versus a position after two time periods of the input excitation () and b time at position x 0.42 m for excitation frequency (i) Hz, (ii) Hz, (iii) Hz and (iv) Hz

A quasistatic homogenization solution (i.e. the simulation with Local Resonance effects ignored) is also shown in order to be able to judge the impact of the Local Resonance coupling. At 200 Hz, the behavior of the LRAM is still captured reasonably well by the quasistatic solution. At 400 Hz, the quasistatic solution cannot reproduce the LRAM behavior at all, showing the pronounced influence of the Local Resonance phenomenon which is picked up remarkably well by the proposed homogenization scheme. Even though the transient simulation is carried out for only 2 time periods without damping, a strong wave decay effect of the stop band is clearly visible in the simulation. At 800 Hz, propagating waves reappear in the LRAM but the wavelength and speed are significantly higher compared to that of the quasistatic wave. This is due to the fact that the resonance of the first eigenmode causes decoupling of the mass from its host matrix reducing the overall effective dynamic mass density of the system. This in turn causes the wave speed to increase resulting in a longer wavelength. Finally, at 1330 Hz, the wave decay phenomena is again active, this time due to the second stop band.

To conclude, the above results illustrate the importance of correctly accounting for the micro-inertia effects in homogenizing the behavior of LRAM. The proposed enriched continuum resulting from the homogeniation and reduction scheme is effective in achieving the correct response.

Conclusions

This paper presented a novel approach towards the modeling and analysis of LRAM with linear elastic constituents. By applying the Craig Bampton reduction to a transient dynamic Computational Homogenization framework, a compact closed form description has been obtained. The resulting equations reveal that an enriched continuum model emerged, with additional degrees of freedom representing the internal dynamics of the microstructure. The method was validated against the DNS solution for a 1D transient compressional wave test of a LRAM structure consisting of a rubber coated lead inclusions at different excitation frequencies. An excellent match was thereby obtained. An order of 3 decade reduction in the problem size with respect to DNS was achieved without any loss of accuracy. The comparison of the results with the quasistatic solution showed the strong influence of Local Resonance on the macroscopic dynamics beyond the low frequency (long wavelength) quasistatic regime.

The main features and advantages of the developed technique, that distinguishes it from other available techniques for the modeling and analysis of materials with Local Resonances (or long wavelength micro-inertia effects in general), can be summarized as follows,In conclusion, the approach presented here aims to pave the path towards integrating core capabilities of Computational Mechanics into the realm of Acoustic Metamaterials. This could provide a breakthrough needed for the practical design and production of such materials for real world applications.

Generality: The full balance of momentum is considered at both scales, which allows the incorporation of complex boundary conditions, transient loading and arbitrary topologies.

Efficiency:

Model reduction allows for a highly compact representation of the microscale problem which naturally implies improved numerical efficiency.

An Enriched macroscopic homogenized continuum description emerges, allowing to use any appropriate solution technique. Therefore discretization schemes with superior convergence properties such as Spectral Element Method etc. can be used as well.

It also offers an additional benefit over standard transient dynamic multiscale implementations in the fact that the homogenization procedure has to only be carried out once and not at each time step.

Extendability: Many possibilities exist for further extending the present formulation, e.g. to incorporate moderate nonlinearities [34, 35], damping, multiphysical effects, material nonlocality (second order gradient) etc.