Computational homogenization of locally resonant acoustic metamaterial panels towards enriched continuum beam/shell structures

Locally resonant acoustic metamaterial (LRAM) panels are excellent candidates for the low-frequency flexural wave ttenuation in thin structures. To enable the efficient analysis and design of LRAM beam/shell structural elements for practical pplications, a computational homogenization method for modelling wave propagation phenomena in LRAM panels is presented n this work. The approach is based on the notion of a relaxed separation of scales, which tailors the methodology to the henomena governed by local resonators embedded in a host medium. The macroscopic LRAM panel is modelled as a thin ontinuum beam/shell described by proper structural kinematics and momentum balance relations. At the microscale, a LRAM nit cell is considered with lamina-wise in-plane boundary conditions and zero out-of-plane tractions, representing the free top nd bottom surfaces of the macroscopic LRAM panel. Under the assumption of linear elasticity, in the relaxed separation of cales regime, the microscale unit cell problem can be represented by the superposition of a static and a dynamic problem, ereby enabling a significant model reduction. As a result, the macroscale effective material properties can be computed once nd for all off-line for a given unit cell. In addition, a new macroscale evolution equation emerges describing the effect of the icroscale internal dynamics on the macroscopic fields through the introduction of new macroscale enrichment variables, which eflect the modal amplitudes of localized microscale modes. The proposed homogenization methodology reveals a high level f accuracy and numerical efficiency compared to the reference direct numerical simulation (DNS) for all relevant analyses: omputation of real and complex dispersion spectra for an infinite LRAM panel, as well as steady-state frequency response nd transient behaviour for a finite LRAM panel. c 2021 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license http://creativecommons.org/licenses/by/4.0/). eywords: Computational homogenization; Locally resonant acoustic metamaterial (LRAM); Beam/shell structural element; Model reduction; nriched continuum; Flexural wave


Introduction
Along with the growing environmental and economic requirements, lighter, but also stronger, designs are being constantly pursued in modern society [1]. Various lightweight designs have emerged with great success. However, most of them are assembled using thin structures, of which the innate high stiffness-to-mass ratio leads to poor performance in terms of noise and vibration isolation [2]. This concern becomes even more pronounced in the low-frequency regime, where in particular the bending modes dominate the vibration behaviour [3,4].
To attenuate waves in thin structures, especially the flexural waves, while preserving their lightweight nature, locally resonant acoustic metamaterial (LRAM) panel designs are being exploited as potential solutions, as sketched in Fig. 1. Different LRAM beams [5][6][7], plates [8][9][10] and shells [11][12][13] have demonstrated capabilities to open wide frequency bandgaps with respect to subwavelength flexural waves, while not requiring large dimensions. Moreover, through adjusting the unit cell material and geometrical parameters, both the lower and upper cutoff frequencies of such stopping bands are tunable, thus providing a high degree of design flexibility and adaptivity [14,15]. The attractive characteristics of LRAM panels are essentially achieved by inducing a localized resonant motion (e.g. originating from the high stiffness contrast between the inclusion and matrix constituents, as investigated in [16] and revisited in [17]) that impedes the out-of-plane linear or rotational motion governing flexural wave propagation, which can be intuitively related to negative effective mass or bending stiffness effects [18][19][20].
Analytical unit cell models can provide remarkable insights for understanding the attenuation mechanisms of flexural waves in LRAM panels [18][19][20]. However, practical LRAM panels may have arbitrarily complex unit cell geometries, macroscopic configurations and boundary conditions. Furthermore, both body (i.e. compressive and shear) and flexural waves may propagate simultaneously and exhibit mutual couplings [11][12][13]. In these cases, analytical models, which rely on many simplifications, are usually inapplicable.
Hence, numerical techniques are frequently employed for fast concept verifications, for which the finite element method (FEM) is one of the most popular approaches. Nevertheless, the computational burden of FEM-based direct numerical simulations (DNS) is heavy in many situations. For instance, a very fine mesh along the thickness direction is typically required to reliably capture bending modes, yielding high computational costs (including the storage cost). Combined with the increasing dimension of the whole macrostructure to be investigated, DNS simulations can rapidly become computationally prohibitive.
Towards addressing the mentioned challenges, homogenization methods may be suitable candidates. Since LRAM panels combine the characteristics of acoustic metamaterials with locally resonant effects and conventional heterogeneous panels for quasi-static applications, representative homogenization methods for LRAMs and heterogeneous panels are separately reviewed in the following.
One notable class of approaches is based on the asymptotic homogenization method [21][22][23]. For the investigation of the wave dispersion relations of an infinite heterogeneous medium, these techniques exhibit comparable accuracy and much higher computational efficiency than a traditional Bloch analysis. In general, the approach revolves around deriving the analytical solutions with respect to a few specific configurations, e.g. a beam lattice assembled by rectangular framed unit cells [21,23] or a layered bi-material containing orthotropic phases [22]. As such, asymptotic homogenization approaches are mathematically rigorous, however, they are usually restricted to simple macroscopic layouts, thereby limiting their application in practice.
These approaches allow for arbitrary configurations and boundary conditions, whereby perfect structural periodicity is not required. In the following, the enriched continuum description-based homogenization method [26] is highlighted. Under the assumptions of linear elasticity and a subwavelength regime, and using a static-dynamic superposition similar to the Craig-Bampton mode synthesis (CBMS) technique [29], the unit cell model can be reduced, providing a closed-form description of a macroscopic continuum enriched with the internal dynamics. Since such a reduction procedure is required only once for a given unit cell design (done off-line), the full multi-scale analysis reduces to a single-scale enriched continuum analysis, yielding a high computational gain. Its accuracy and effectiveness in dispersion analyses, frequency domain analyses and transient analyses have been demonstrated in [26,30,31]. It must be pointed out that these dynamic computational homogenization frameworks are developed using a classical continuum solid description at both micro and macroscales. As a result, they may not be directly applied to study thin LRAM panels with the (mostly) single-layer configurations, (approximately) free top and bottom surfaces, and (possibly) simultaneous coupled propagation of both body and flexural waves.
Limited attempts to homogenize panels with internal dynamics have been carried out, able to describe the interaction with the subwavelength flexural wave propagation. To this end, the effective medium theory [32][33][34] and asymptotic expansion [35][36][37][38] are often employed, yielding analytical homogenization solutions for certain designs or conditions. The applicability of some of the approaches [33,34] is limited to homogeneous panels with periodically attached local mass-spring resonators. Periodic reticulated beams made of symmetric unbraced framed cells were investigated in [35], bi-laminate periodic beams were considered in [36], and the derivations in [37,38] were established for unidirectional or bidirectional periodic ribbed plates with local resonance. The specific formulations to investigate near and far-field scattering phenomena were discussed in [32], where the localized modal terms need to be known in advance. Recently, a homogenization method based on the diffusion matrix model has been developed in [39], where the dynamic effective material properties are determined numerically by minimizing the reflection coefficients at the interface between the unit cell and equivalent model, thus allowing for a complex geometry. Nevertheless, the accuracy of this homogenization method largely relies on the proper choice of modes, which are often unknown under general loading conditions. To the authors' best knowledge, no generalized continuum description based homogenization framework for LRAM panels has so far been presented.
As stated early, the successful development of computational homogenization frameworks for LRAM solids has benefited from extensions of the traditional computational homogenization frameworks for heterogeneous solids under quasi-static conditions. A natural point of departure to homogenize LRAM panels are methods developed for quasi-static heterogeneous panels [40][41][42][43]. An early general multi-scale framework can be found in [40], applicable to thin and thick sheets, arbitrary micro and macroscopic layouts, and boundary conditions. After that, a variety of improvements and modifications have been proposed, e.g. those compatible with standard beam/shell FE formulations [41], not requiring high-cost scale-coupled computations [43] and specifically reduced for layered or stacked composite panels [42].
Among the above computational homogenization frameworks, the one presented in [41], which is established under the macroscopic continuum beam/shell and microscopic continuum solid descriptions is selected as the point of departure in this work to develop a general computational framework for LRAM panels. Here, it will first be extended to the dynamic regime by taking into account the inertial effects at both scales and consistently reformulating the scale transition relations. Next, using the assumptions of linear elasticity and a subwavelength regime, a static-dynamic superposition [26] will be adopted, leading to a closed-form description of the homogenized continuum beam/shell at the macroscale, enriched by micro-inertia. The computational homogenization framework proposed in this work (i) is applicable for arbitrarily complex micro-and macroscopic geometries and boundary conditions and (ii) can capture both body and flexural waves in the subwavelength regime accurately, all at a low computational cost. This paper is organized as follows. The continuum beam/shell description is established in Section 2, which is adopted later as the macroscopic description of the heterogeneous panel under dynamic loading. Section 3 summarizes the microscale problem description. In Section 4, the separation of scales conditions is specified that define the applicability range of the proposed approach. Combining the introduced key concepts, the dynamic computational homogenization framework for LRAM panels is formulated in Section 5. For the linear elastic case, in Section 6, the discretized unit cell model is reduced using a static-dynamic decomposition and mode synthesis. Combining this reduced unit cell model with the upscaling relations, finally, yields the enriched macroscopic continuum beam/shell description. In Section 7, the dispersion spectra of an infinite LRAM beam are obtained from the homogenization-based dispersion analysis and compared to those obtained by a direct application of the Bloch analysis. In Section 8, the steady-state and transient macroscopic responses of a finite LRAM beam are computed by solving the enriched macroscale problem, which is compared to the full-scale FEM-based DNS. The computational cost reduction provided by the homogenization and model reduction for different analysis types are discussed as well. Finally, the conclusions are given in Section 9.
Scalars, vectors, 2nd-order tensors, 3rd-order tensors and 4th-order tensors employed throughout this paper are denoted as e.g. a, ⃗ a, A, A (3) and A, respectively. The following notation for vector and tensor product operations are employed together with Einstein's summation convention: the dyadic product ⃗ a ⃗ b = a i a j ⃗ e i ⃗ e j , the dot product A · B = A i j B jk ⃗ e i ⃗ e k , and the double dot product A : B = A i j B ji , with {⃗ e 1 , ⃗ e 2 , ⃗ e 3 } the Cartesian vector basis. The column and matrix assemblies are, denoted by (•) and (•), respectively. The inverse and transpose of a 2nd-order tensor or a matrix are denoted by (•) −1 and (•) T , respectively. The left transpose of a higher-order tensor is denoted by (•) LT . The macro and microscopic quantities are distinguished by the subscripts "M" and "m", respectively.

Macroscale: continuum beam/shell description
In this section, a general description of the continuum shell used at the macroscale is given. 1 This shell description is rather standard and can be found in many text books (see e.g. [44]); it is summarized here for the sake of completeness and in order to define a proper setting for the subsequent homogenization framework. The continuum shell kinematics is first introduced, followed by the definition of the through-thickness stress resultants and the weak and strong forms of the balances of the linear and rotational momenta. In this section, the subscript "M" is temporarily omitted since only one (macro) scale is relevant here.

Kinematics
Fig. 2 schematically illustrates the initial ( Fig. 2(a)), current ( Fig. 2(c)) and rigid rotation neutralized ( Fig. 2(b)) configurations of a macroscopic continuum shell. Here, for simplicity, the initial configuration is depicted flat, but the subsequent formulation allows to account for an initial curvature as well. In order to quantify the kinematics, a reference surface is defined in the middle between the top and bottom surfaces. All material points with the same distance from the reference surface constitute a lamina. Each material line in the initial configuration normal to the midsurface is defined as a fibre. In addition, a co-rotational laminar basis {⃗ e c 1 ,⃗ e c 2 ,⃗ e c 3 } is adopted (marked by the superscript "c"), such that the (⃗ e c 1 ,⃗ e c 2 )-plane is tangent to the reference surface and⃗ e c 3 , in the following simply denoted as⃗ e c , is normal to this tangent plane. The co-rotational laminar basis rotates rigidly with the continuum shell frame and can thus be regarded as a local Cartesian vector basis. With respect to the co-rotational laminar basis, the out-of-plane and in-plane components of a vector are represented by (⃗ •) and (⃗ •), respectively, such that (⃗ •) = (⃗ •) ·⃗ e c⃗ e c and (⃗ •) = (⃗ •) · (I −⃗ e c⃗ e c ), where I is the 2nd-order identity tensor. This notation is also extended to tensors, constructed through the dyadic product of the in-plane or out-of-plane projected vectors, e.g.
For convenience, the subsequent derivations will be performed between the initial and rigid rotation neutralized configurations. For problems involving large rigid rotations, corresponding expressions can be readily obtained using the standard pull-back relations. Position vectors of a material point in the initial and rigid rotation neutralized configurations are written, respectively, as where the subscript "r" denotes the reference surface; ⃗ D is the director in the initial configuration, which is an upward directional vector of fibre, coincident with the normal ⃗ N to the reference surface; ⃗ d denotes the director in the rotation neutralized configuration, which is, in general, non-coincident with the normal ⃗ n, thus enabling description of fibre length change and transverse shear effects; η is the out-of-plane coordinate along the fibre, belonging to the thickness domain H = [−H/2, H/2], with H the total shell thickness.
Subtracting Eq. (1) from Eq. (2) gives the displacement field as where ⃗ u r reflects the reference point displacement induced by the linear motion of the reference surface and η ⃗ b the bending displacement induced by the fibre rotation, with Considering thin structures only, the subsequent formulation will be restricted to thin continuum shells, whose thickness is much smaller than the out-of-plane radius of curvature. Therefore, the Kirchhoff-Love theory (corresponding to the Euler-Bernoulli theory for continuum beams) [45] is applicable here, i.e. all fibres are assumed to always stay straight, inextensible and normal to the midsurface during the deformation process. As a result, for a thin shell, the following kinematic relation exists between the bending-related rotation vector ⃗ b and the out-of-plane displacement of the reference surface⃗ u r : where⃗ ∇ 0 = ∂(•)/∂ ⃗ X r is the gradient operator with respect to ⃗ X r . The deformation of an infinitesimal material line element between the rigid rotation neutralized and initial configurations can be written as with the deformation gradient tensor following from Eqs. (1) and (2) as Here is the gradient operator with respect to ⃗ X . Using the previously defined decomposition, the in-plane and out-of-plane components of the infinitesimal material line element in Eq. (5) can be written as (see [41] for more details) which are the in-plane deformation gradient tensor, deformation curvature tensor and in-plane gradient of transverse position component, respectively. Note, that in the rigid rotation neutralized configurationF =0 and therefore will not be included in the following derivations. The transverse shear deformation termF =⃗ d⃗ e c has been neglected, according to the thin shell assumption that each fibre always stays normal to the reference surface.

Stress and momentum resultants
The kinetics of a shell continuum is typically described by the through-thickness resultant quantities. Denoting the 1st Piola-Kirchhoff stress tensor by σ and the linear momentum rate by⃗ p, the resultants are given by which denote the stress resultant, linear momentum rate resultant, couple-stress (moment) resultant and rotational momentum rate resultant, respectively.

Weak and strong forms of the balance equations
The weak form for a thin continuum shell, i.e. neglecting the transverse shear contributions, can be stated as follows: find (⃗ u r , ⃗ b) such that for all admissible test functions (δ ⃗ u r , δ ⃗ b), the balance holds between the internal, kinetic and external virtual powers: Here, A r denotes the domain of the reference lamina (which is a surface) in the initial configuration; C r is the external boundary of A r , with the unit outward normal ⃗ N A and its in-plane component⃗ N A . Using kinematic relation (4) between⃗ b and⃗ u r , applying the chain rule and Gauss divergence theorem and taking into account that the weak form should hold for all admissible δ⃗ u r and δ⃗ u r , leads to the strong form of the momentum balance equations for a thin shell: with the traction and couple-traction resultant boundary conditions on the external contour C r given bŷ The macroscopic balance equations need to be complemented by constitutive relations connecting the resultantform responsesN,M,⃗ P and⃗ Q to the kinematic quantities. In this work, these relations will not be pre-assumed for the macroscopic continuum shell. Instead, they are to be obtained through consistent scale transitions, following the framework of the dynamic computational homogenization, which will be elaborated later. As a final remark, note that the resultant momenta⃗ P and⃗ Q are treated here as a constitutive quantities, instead of directly expressing them in terms of the (effective) mass density and accelerations, which cannot be done trivially for metamaterials.

Microscale: continuum solid description
To each macroscopic fibre, a microscopic unit cell is associated, as sketched in Fig. 3. Here, the unit cell thickness H spans the whole thickness of the macroscopic beam/shell, while the unit cell in-plane dimensions L m and W m are representative for the considered microstructure (e.g. a repeating unit). Notice, that although the thin beam/shell theory employed at the macroscale requires H to be much smaller than the in-plane dimensions of the whole macrostructure, this does not apply to the unit cell, where H is typically of the same order as L m and W m . In addition, it is assumed that the unit cell out-of-plane radius of curvature is much larger than its in-plane dimensions, both in the initial and deformed configurations. Hence, as along as this condition holds true, a small out-of-plane curvature of the unit cell is allowed. Based on this assumption, the co-rotational laminar basis {⃗ e c 1 ,⃗ e c 2 ,⃗ e c 3 } at the corresponding macroscopic point can be used as a Cartesian vector basis for the respective microscopic unit cell.
Let ⃗ X m and ⃗ x m denote, respectively, the initial and deformed position vectors of a material point in the unit cell domain V m . The deformation within the unit cell is described by with the deformation gradient tensor The balance of linear momentum, disregarding the body forces, is expressed as with σ m the stress and ⃗ p m = ρ m⃗ u m the momentum. The corresponding weak form that holds for all admissible test where S m is the external surface of V m , with the outward normal indicated by ⃗ N m . The constitutive stress-strain relations for all microstructural constituents are assumed to be known. The dynamic computational homogenization framework to be established here allows for general nonlinear and time-dependent material laws, in a general form written for each microscale constituent (α) as where t is the time. The boundary conditions required to complete the microscale boundary value problem will be determined through the macro-to-micro scale transition, i.e. downscaling, which is elaborated in Section 5.

Separation of scales
The recently developed transient computational homogenization methods [24,26,30,31,46] for 3D metamaterials operating in the subwavelength regime rely on a relaxed separation of scales, where a wavelength contrast between the inclusion and matrix constituents is introduced (i.e. accounting for both the mass density contrast and the elastic stiffness contrast, see e.g. [17]). In particular, the matrix host medium is required to be in its long wavelength regime, while the wavelength in the resonating inclusions is allowed to be comparable to the inclusion size, thereby allowing for non-negligible micro-inertia effects. Note, that the inclusions within LRAM are typically designed with low stiffness and heavy mass, and placed sufficiently far from each other, to avoid mutual interaction and thereby trigger the local resonance.
The above relaxed scale separation definition is now specialized for LRAM panels, where the flexural wave propagation is of particular concern. Note, that the transverse body waves are negligible due to the small thickness. Accordingly, a lamina-level relaxed separation of scales can be defined as where the subscripts "mat" and "inc" denote the matrix and inclusion, respectively; the superscript "A" denotes the lamina; h and l are the representative transverse thickness and in-plane length of a constituent, respectively; λ bo is the characteristic body wavelength while λ fl the characteristic flexural wavelength within a constituent or the whole macrostructure, under the considered excitation. The requirements on l A in each lamina ensure that the considered body and flexural waves are in their subwavelength regimes, which enables the homogenization. The additional requirements on h ensure the applicability of the thin beam/shell theory. The lamina-level separation of scales in Eq. (18) defines the range of applicability of the dynamic computational homogenization framework for LRAM panels, to be presented next.

Dynamic computational homogenization framework for LRAM panels
Based on the established macroscopic beam/shell and microscopic continuum solid descriptions, and the definition of the lamina-level separation of scales, the dynamic computational homogenization framework for LRAM panels is formulated. The generalized macroscopic deformation mapping is imposed on the unit cell through the macro-to-micro scale transition, resulting in appropriate lamina-wise boundary conditions. Upon solution of the microscale boundary value problem, the macroscopic resultant-form responses are obtained through the micro-to-macro scale transition.

Macro-to-micro scale transition
In the following, the deformed configuration of the microscopic unit cell corresponding to the rigid rotation neutralized macroscopic shell configuration will be considered. Applying the macroscopic in-plane deformation gradient tensorsF M andK M to the microscopic unit cell in accordance with Eq. (7), with account thatF M =0, yields the lamina-wise microscopic relative position vector field, given in terms of its in-plane and out-of-plane components by where ⃗ X m,r denotes the initial configuration position vector of a reference point (marked by the subscript "r"), which is located within the reference (i.e. middle) lamina of the unit cell to match the kinematic description of the macroscopic fibre; the position vector of this point after deformation is denoted by ⃗ x m,r ; ⃗ w m reflects the microfluctuation field induced by the heterogeneities. The microscopic displacement field of the unit cell can be represented by . In the following, the origin of the microscopic Cartesian vector basis will be introduced at ⃗ X m,r and taken to coincide with the associated macroscopic co-rotational laminar basis.
In the computational homogenization frameworks for continuum solids [47,48], the unit cell volume average of the microscopic deformation gradient tensor F m is equated to the macroscopic deformation gradient tensor F M . In shell homogenization [41], the unit cell is associated to a through-thickness material line (fibre), thus requiring F A M (η) along the macroscopic fibre to be equal to the microscopic lamina average of F A m (η). This homogenization requirement yields the downscaling relations for a continuum shell: Here C m denotes the closed contour of A m with the unit outward normal ⃗ N m and its in-plane component⃗ N m . Without loss of generality, if ⃗ X m,r is set at the unit cell centre and assuming a small out-of-plane curvature compared to the unit cell dimensions, the last term in Eq. (23b) vanishes Finally, with account for Eq. (24), Eq. (23) gives the lamina-wise constraints on ∆ ⃗ w m : Constraints (25) are usually too weak and can be satisfied using several stronger assumptions, e.g. prescribing zero or periodic microfluctuations on each lamina boundary, resulting in the uniform displacement or periodic boundary conditions, respectively. It has been extensively confirmed in the literature, e.g. [49], that the periodic boundary conditions provide a better convergence of the effective quantities versus the unit cell size. Hence, they are also adopted here. Note, that the periodic boundary conditions are violated in general for the dynamic cases. Instead, the Bloch-Floquet boundary conditions hold for the wave propagation in an infinite medium. Nevertheless, the periodic boundary conditions, that correspond to Bloch-Floquet conditions in a long wavelength regime, are acceptable here, based on the separation of scales defined in Eq. (18), which requires the long wavelength condition for the matrix containing the inclusions.
Since the unit cell has an in-plane geometrical periodicity on each lamina boundary, the in-plane normal vectors are anti-periodic. Then, Eq. (25) is automatically satisfied when lamina-wise periodicity of the micro fluctuation field is imposed: where ξ specifies the in-plane local coordinate along the lamina contour C m and the superscripts "−/+" denote opposite transverse boundary pair. Note, that enforcing periodic boundary conditions on the microfluctuation field leads to anti-periodicity of the lamina-wise tractions, which follows from the condition of zero work by the constraints. This property will be exploited in the upcoming micro-to-macro scale transition. Using Eq. (19) in combination with Eq. (20), constraints (26) can be expressed in terms of displacements on the microscopic lamina boundaries (i.e. the transverse surfaces) as functions of the prescribed macroscopic kinematic quantities ⃗ u M,r , ⃗ b M ,F M andK M . Since a plane stress state is assumed for the macroscopic beam/shell in the thickness direction, zero traction boundary conditions are prescribed on the top and bottom boundaries of the unit cell (marked by the superscripts "U" and "D", respectively): This completes the microscale boundary value problem formulation.

Micro-to-macro scale transition
After solving the microscale boundary value problem, the microscopic stress and momentum can be readily homogenized towards the macroscale, based on the Hill-Mandel principle [50]. The classical Hill-Mandel condition for a continuum solid requires the volume average of the unit cell virtual work to be equal to the virtual work per unit of volume at the associated macroscopic material point. For a macroscopic beam/shell, the stress, couple-stress and momentum are described in their resultant forms. Therefore, the Hill-Mandel condition needs to be formulated in a compatible form [41], i.e. by equating the unit cell virtual power per unit reference lamina area, δW A m , to the virtual power per unit of area at the associated macroscopic fibre δW A M : The left-hand side of Eq. (28) is given by with A m,r = A m (η = 0). To obtain the second equality in Eq. (29), the chain rule with account for the microscopic balance of linear momentum Eq. (15) and Gauss's theorem have been used.
In accordance with Eq. (20), δ ⃗ u m can be expressed as with δ∆⃗ x m = δ∆⃗ x m + δ∆⃗ x m , revealing the in-plane and out-of-plane components that follow from Eq. (19) as whereS m denotes the transverse boundary of the unit cell. The macroscopic virtual work per unit of volume area, i.e. the right-hand side of Eq. (28) naturally follows from Eq. (10) and reads Equating (32) and (33) leads to the upscaling relations for the macroscopic resultants: Eq. (34) enables a straightforward calculation of the macroscopic resultant responses, using quantities from the unit cell transverse boundaries only. Note, that different from the quasi-static problem (where the exact position of the reference point ⃗ X m,r is irrelevant due to the equilibrium of the boundary tractions), the macroscopic stress state in the dynamic problem is dependent on the choice of ⃗ X m,r . Moreover, both⃗ t m and⃗ t m onS m , in general, contribute toM M . The transverse traction-related term in Eq. (34b) essentially reflects the relative sliding between the laminae induced by the transverse shear effects at the scale of the unit cell. However, for the problems where the transverse shear deformation effects can be neglected (which is the case here), this term can be disregarded, implying thatM M is only determined by⃗ t m . With this simplification, by applying the lamina-wise Gauss's theorem, the chain rule and Eq. (15), the macroscopic resultants (34) can be expressed in terms volume integrals, readinĝ Note the existence of a coupling between the macroscopic stress resultants and microscopic inertial effects in Eq. (35). The right hand side of Eq. (35) also illustrates that the macroscopic resultant responses can be obtained through the homogenization of the microscopic responses in each lamina, followed by the integration through the thickness, indeed consistent with the definition of the macroscopic resultants, as introduced in Eq. (9).

Model reduction towards an enriched macroscopic continuum beam/shell
The transient computational homogenization for the dynamic heterogeneous panels presented in the previous sections can be employed in general cases, including nonlinearities. If small deformations and linear elasticity are assumed, the coupled two-scale framework can be further reduced to a single-scale continuum beam/shell framework, enriched with extra micro-inertia field variables. This model reduction is elaborated in this section. First, the microscopic balance equations of linear momentum are discretized in space, followed by the application of the lamina-wise periodic boundary conditions. Next, the static-dynamic decomposition and mode synthesis are adopted for further reduction of the unit cell model. Finally, the application of the upscaling relations to the reduced model results in the macroscopic governing equations, describing an enriched single-scale continuum beam/shell.

Finite element discretization of the microscale problem
The established microscale boundary value problem can be solved using any appropriate method. Here, a FEMbased implementation is used. In the following, the discussion will be restricted to small-strain linear elastic behaviour for all microscopic constituents, which is an assumption that is usually valid for LRAM panel designs and analyses. Therefore, for the subsequent derivations the microscale constitutive laws Eq. (17) are specified for each microstructural constituent (α) as with C m the 4th-order elasticity tensor, which can be anisotropic, and ε m the infinitesimal linear strain tensor.
Applying the standard FE discretization to Eq. (16) yields the discretized microscale balance equations of linear momentum, reading where B m and W m denote the assembled microscopic stiffness and mass tensor matrices, respectively; ⃗ ũ m and ⃗ f m are the nodal displacement and force vector columns, respectively.
Eq. (37) has to be complemented with boundary conditions, namely the zero-traction boundary conditions on the top and bottom surfaces, Eq. (27), the lamina-wise periodic boundary conditions, Eq. (26), and the conditions prescribing the translation and bending related rotation of the macroscopic fibre. The latter are imposed on the microscopic unit cell through the prescribed displacements of several "prescribed" nodes (marked by the subscript "p"), that from Eqs. (19) and (20), can be written in a general form aŝ where η m,p is a diagonal matrix containing the through-thickness coordinates of the prescribed nodes η m,i ; ∆⃗ X m,p is the column with the (relative) positions of the prescribed nodes ∆⃗ X m,i and ∆X m,p is the column with the dyadic products of ∆⃗ X m,i , i.e. its ith element corresponding to the ith prescribed node is defined as ∆X m,i = ∆⃗ X m,i ∆⃗ X m,i ; andĤ M =F M −Î. The particular selection of the prescribed node set for a unit cell that represents for a macroscopic continuum beam is discussed in Appendix A.
The lamina-wise periodic boundary conditions, Eq. (26) can be rewritten in terms of ⃗ u m instead of ∆⃗ w m , as detailed in Appendix A, making them compatible with a standard FEM framework. In a general form, they can be implemented as a set of linear constraint relations between the displacements of the (retained) "independent" and (tied) "dependent" nodes (marked by the subscripts "i" and "d", respectively), reading as where T m,di denotes the dependency matrix containing the corresponding coefficients.
Partitioning the system of Eqs. (37) into the dependent and independent degrees of freedom (DOFs) sets, and eliminating the dependent DOFs from Eq. (37) using Eq. (39), leads to the system of equations written in terms of the independent DOFs only with Here the superscript "*" denotes the quantities after elimination of the dependent DOFs.

Static-dynamic decomposition and model reduction
In accordance with the relaxed separation of scales (see Section 4), the matrix behaviour is in a long wavelength regime (i.e. it can be approximated as quasi-static) and thereby the unit cell boundaries can be regarded as sufficiently stiff. Consequently, the Craig-Bampton mode synthesis (CBMS) technique [29] can be adopted to obtain a reduced description of the total dynamic response of the unit cell under the prescribed boundary kinematics [26].
The reduced model is constructed using two reduced bases: the quasi-static modes and fixed-boundary dynamic modes, as schematically shown in Fig. 4, where a 2D unit cell (i.e. corresponding to a continuum beam) is depicted for illustration purposes only; the subsequent derivations remain valid for the general 3D case. The quasi-static deformation modes are computed by neglecting the inertial effects. The fixed-boundary dynamic modes are a set of internal dynamic eigenmodes for which the prescribed nodes are fixed; the dynamic modes contain the information related to the localized resonance within the unit cell. In the following, the static and dynamic modes will be elaborated separately and then superposed to obtain the total solution.
For the subsequent derivations, the independent nodes are further subdivided into the "prescribed" and "free" node sets (marked by the subscripts "p" and "f", respectively). Hence, Eq. (40) can be partitioned as where it has been taken into account that ⃗ f * m,f = ⃗ 0, since no external force exists on the free nodes after imposing the tying constraints.

Quasi-static deformation
After omitting inertial terms in Eq. (42), its solution can be expressed in terms of the displacements of the prescribed nodes as with the static condensed tensor matrix (also known as the Schur complement) given by In Eq. (43), the superscript "qs" denotes the quasi-static deformation and I m,pp is the identity tensor matrix, such that ⃗ ũ qs m,p = I m,pp · ⃗ ũ m,p .

Internal dynamics
The internal dynamics modes account for the inertial effects were disregarded above for the quasi-static subproblem. The solution for the internal dynamics subproblem is obtained as a superposition of eigenmodes, computed by fixing the prescribed nodes, imposing the zero traction boundary conditions on the top and bottom boundaries and the lamina-wise periodic boundary conditions on the transverse boundaries. Note, that since the matrix is in the long wavelength regime, the Bloch-Floquet boundary conditions would degenerate to the adopted periodic boundary conditions. The solution in the dynamic regime can be represented as where the superscript "dy" denotes the internal dynamics, ⃗ 0 is the matrix of zero vectors, enforcing ⃗ ũ dy m,p = ⃗ 0; ⃗ φ is a matrix, each column of which contains the normalized modal displacement vector field ⃗ φ n , n = 1, 2, 3, . . . , N d , with N d the number of modes included in the analysis; ζ is the column containing the corresponding modal displacement amplitudes ζ n , which are not known at this stage and will only be computed at the macroscale after the two-scale coupling. To obtain ⃗ φ n , harmonic solutions are sought for Eq. (42), leads to a series of eigenvalue problems with the normalization condition given by In Eq. (46), ω n denotes the eigenfrequency (sorted in an ascending order) corresponding to ⃗ φ n . Note, that in the frequency range of interest, it is usually unnecessary to account for all eigenmodes and only a few relevant locally resonant modes are important for constructing the reduced basis. A simple eigenmode selection criterion will be introduced later.

Linear superposition
Additively superposing Eqs. (43) and (45) reconstructs the total unit cell response as which reveals that the original solution for ⃗ ũ m,f is now reduced to the solution to that for ζ, only containing a very limited number of unknowns (N d ).
Substituting Eq. (48) into Eq. (42), 2 and taking into account Eqs. (46) and (47) yield a reduced dynamic problem as with In Eq. (49), the superscript "uζ " in particular denotes the coupling between the internal dynamics and the total unit cell response; ω 2 n is a diagonal matrix containing ω 2 n .
After substituting Eqs. (38) and (49a) into Eq. (51) and using the orthogonality between the in-plane and outof-plane components, the homogenized constitutive relations for the reduced unit cell model can be identified, as the following set of compact closed-form expressions (the detailed derivation can be found in Appendix B): 2 For this operation, the general transformation rule of a column of quantities ũ to ũ * , related by the transformation matrix T as ũ = T ũ * has been applied to a linear system of equations K ũ = f, leading to a transformed system of equations K * ũ * = f * with K * = T T K T and f * = T T f.
with the constitutive quantities given by Here, C N M and C M M are the quasi-static effective membrane and bending stiffness tensors; D N M and D M M are the effective membrane and bending elastic inertia tensors; ρ P M and ρ Q M are the generalized quasi-static effective mass density tensors; the tensorial and vectorial coefficients h N n , h M n , ⃗ j P n and ⃗ j Q n reflect the coupling between nth modal DOF ζ n and the macroscopic stress, couple-stress and momentum rate resultants. Each mode introduces at least one non-zero ⃗ j n or h n (which is symmetric as long as no locally rotational coupling takes place). Note, that each localized mode can, in general, contribute to both body and flexural waves simultaneously.
Finally, it should be emphasized that, as opposed to many other homogenization methods for (locally resonant) metamaterials, all effective material properties obtained here are frequency independent. This gives a significant advantage when solving problems in the time domain.

Selection of dynamic modes
The extracted effective material properties can be used to define a criterion for the selection of (a minimal number of) eigenmodes for the reduced dynamics basis. To this aim, first, the 2nd-order modal mass density tensor J n and the 4th-order modal inertia density tensor H n are defined as J n = ⃗ j n ⃗ j n /A m,r and H n = h n h n /A m,r , which allows to define a generalized modal mass and elastic inertia fraction for each mode n following [30], given by where N dim is the dimension index. To achieve an adequate accuracy, the sum of all modal fractions of the localized modes should approximately be equal to the total generalized mass and elastic inertia fractions of the inclusion, µ inc , which can be evaluated beforehand by where ρ inc and D inc can be computed similarly to ρ M and D M using Eqs. (53c)-(53f), by setting the contribution of mass matrices for matrix material to zero when assembling corresponding W uu m . Thus, a mode selection criterion can be established as where ϵ ρ tol and ϵ D tol are predefined tolerances. An example of this mode selection procedure will be given in Section 7.

Governing equations for the enrichment field variables
As mentioned early, the enrichment field variables ζ n need to be solved for at the macroscale. Substituting Eq. (38) into Eq. (49b) upscales the balance equations of the internal dynamics, leading to (the derivation can be found in Appendix B) where h N n , h M n , ⃗ j P n and ⃗ j Q n are given in Eqs. (53g)-(53j). Eq. (58) is an ordinary differential equation (ODE) for each dynamic mode, to be solved at the macroscale coupled with the momentum balance Eqs. (11), which are partial differential equations (PDE).

Summary of the reduced homogenization framework
To summarize, the developed reduced computational homogenization framework leading to a homogenized enriched continuum at the macroscale consists of the following steps. First, during the off-line stage, the unit cell is discretized (e.g. by using finite elements). Upon assembly, the static condensation and eigenfrequency analyses are performed and the frequency independent constitutive tensors are computed according to Eq. (53). This off-line stage has to be performed only once for a given unit cell. The on-line analysis then consists of solving the macroscale boundary value problem Eqs. (11)- (12), with the constitutive relations (52), coupled with the internal dynamics evolution Eq. (58). As a result, the fully coupled two-scale computational homogenization analysis reduces to an enriched single-scale analysis with a limited number of enrichment variables ζ n corresponding to the number of relevant locally resonant modes in the frequency interval of interest. This results in a significant reduction of the computational cost. A schematic overview of the enriched homogenization framework is shown in Fig. 5.

Numerical simulations on an infinite LRAM panel
The developed computational homogenization method is next applied to the analysis of the steady-state wave propagation in an infinite LRAM panel with negative effective mass effects. Using the plane wave transformation, the dispersion eigenvalue problem for the enriched continuum beam/shell is obtained and the accuracy of the computed dispersion spectrum is compared to the reference solution from Bloch analysis. A procedure for a priori establishing frequency limits of the applicability of the homogenization method is developed and tested.
In the examples considered in the following, the rotational inertia effects are negligible, and, therefore, without loss of generality, the terms related to⃗ Q were disregarded in the implementation.

Dispersion relation
The plane wave transformation in an infinite medium reads where ⟨•⟩ denotes the transformed quantity, i is the imaginary unit, ω the angular frequency, ⃗ k = k ⃗ e w the wave vector, k = ∥ ⃗ k∥ and ⃗ e w the wave direction. Inserting the plane wave expression Eq. (59) into Eqs. (11), (52) and (58), taking into account thatĤ T M =⃗ ∇ 0⃗ u M,r andK T M =⃗ ∇ 0⃗ b M = −⃗ ∇ 0⃗ ∇ 0 (⃗ e c ·⃗ u M,r ), leads to the eigenvalue problem: with whereê w =⃗ e w⃗ e w . In this formulation, the wave vector ⃗ k is a parameter. Solving Eq. (60) for angular frequencies ω provides the dispersion relation in the so-called ω(k) form, which reveals only the real branches of the dispersion spectrum. From Eqs. (60)-(61) it is also clear that only in-plane propagating waves can be supported, as expected for a thin beam/shell.
To find the evanescent wave branches as well, the k(ω) form of the eigenvalue problem is further derived by eliminating ⟨ζ⟩ from Eq. (60), resulting in: with in which the generalized dynamic effective mass density and elasticity tensors (marked by the subscript " ") are given bŷ with J and H the tensors already defined in Eq. (54). The 3rd-order dynamic cross-coupling coefficient tensor for the flexural wave R (3) is defined as The above eigenvalue problem in the k(ω) form can also be written as two uncoupled subproblems for in-plane and out-of-plane modes, respectively: which demonstrates that the dynamic effective mass densityρ P and elasticity tensor C N relate to the body wave, while the flexural waves are governed byρ P and C M , as well as R (3) .

Case study of an infinite LRAM panel
The real and complex dispersion spectra of the homogenized enriched continuum computed according to Eqs. (60) and (62), are compared to the reference solution provided by the Bloch analysis on the fully discretized unit cell. The Bloch analysis specialized to a unit cell representing a continuum beam or shell, i.e. having free top and bottom surfaces was discussed before in [13]. To obtain the k(ω) form, a procedure similar to the one presented in [51] is followed. To facilitate the comparison of the computational costs, both the reference Bloch analysis and the homogenization-based dispersion analysis have been implemented using MATLAB, on a standard desktop computer with an Intel (R) Core (TM) i5-4200H CPU, 2.80-3.40GHz Clock Speed and 4GB RAM. No parallel computation is used.

Geometric model configurations
At the macroscale, an infinitely long LRAM panel with free top and bottom surfaces is considered, see Fig. 6(a). A plane strain condition is assumed in the⃗ e c 2 -direction; in this case, the continuum beam description can be used for the homogenized panel. The unit cell follows a classical negative mass design [52] consisting of a relatively stiff, but light matrix, a soft coating and a heavy core, see Fig. 6(b). The material properties of the unit cell constituents and the geometrical parameters used in the simulations are given in Tables 1 and 2, respectively. The unit cell is discretized with 2D plane strain linear quadrilateral elements, with an element size in the range 0.1-0.5 mm, leading to a total number of 1554 DOFs. This finite element mesh has been confirmed to yield the converged solution for the eigenfrequency analysis in the considered (low) frequency range.

Off-line stage: effective material properties
The constitutive quantities of the considered unit cell have been computed according to Eq. (53) and are given in Table 3. The first six modes with the lowest eigenfrequency, shown in Fig. 7, have been included into the dynamic reduced basis. The total in-plane and out-of-plane modal mass fractions of these six modes are evaluated as 44.4% and 44.2%, respectively, which have been considered to be sufficiently close to the inclusion mass fraction of 46.6%. The total modal linear elastic inertia fraction is negligible, in agreement with the small linear elastic inertia fraction of inclusions 5.17%; the total modal rotational elastic inertia fraction is 54.6%, well matching the inclusion rotational elastic inertia fraction 57.0%. Based on these results and following the eigenmode selection procedure, Eq. (57) with a choice of ϵ ρ tol = 5%, N d = 6 has been selected for the subsequent analysis. Note, for the considered unit cell design, where negative mass effects are dominating, the elastic inertia plays a minor role and thus the choice of ϵ D tol is not important.

On-line stage: dispersion analysis
The real dispersion spectrum of the considered infinitely long LRAM panel computed using the ω(k) approach, Eq. (60), is shown in Fig. 8(a) together with the reference solution, accompanied by the modal shapes from the Bloch analysis at the normalized wave number kλ ref /2π = 1, with λ ref = 5L m (≈ 6H ). As expected, the localized modes indeed correspond to the modes used for the construction of the reduced unit cell dynamic basis (see Fig. 7). Fig. 8(b) shows the components of the generalized dynamic effective mass density tensor ρ P , computed in accordance with Eqs. (64a) and (64b). Fig. 8(a) indicates that both flexural and body wave frequency bandgaps emerge in agreement with the corresponding negative effective mass intervals, resulting from the out-of-plane and in-plane localized resonances, respectively. In the frequency ranges where these two bandgaps overlap, the full frequency bandgaps are opened on both body and flexural wave branches. In addition, it can be noticed that two localized rotational modes (1 and 4) do not contribute the bandgaps, since, for the thin LRAM panel considered here, the unit cell rotational modes do not couple to the body and flexural waves. 3 Comparing the real dispersion spectra obtained using both approaches reveals that in the low-frequency regime considered (in accordance with the separation of scales conditions Eq. (18)), the homogenization-based dispersion analysis excellently matches the reference Bloch analysis. The maximal deviation on the phase velocity equals 1.61%.  Next, the evanescent wave properties are investigated. Note, that the total number of eigenmodes to be solved in the k(ω) form of the reference Bloch analysis and the homogenization-based dispersion analysis are different, since this number is determined by the order of the eigenvalue problem and number of unknowns, both of which are different for the two methods. For the reference Bloch analysis, the problem is 2nd-order and, for the 2D unit cell discretization shown in Fig. 6(b), involves 22 unknowns, resulting in 44 solutions, while for the homogenizationbased analysis, Eq. (62), the problem is 4th-order and has 2 unknowns (for the beam problem considered here), leading to 8 solutions for each frequency. Nevertheless, only the relevant branches (two for the body wave and two for the flexural wave) are of interest here. These have been used to plot the complex dispersion and phase angle spectra of the considered infinitely long LRAM panel in Fig. 9.
The powerful capability of the homogenization method is validated on capturing not only the propagative but also evanescent waves, even in the frequency range where no negative effects occur. The well-known phenomenon of the appearance of evanescent flexural wave branches in thin structures [4] is well picked up by the homogenization method, as demonstrated in Fig. 9(a). A small difference is visible in Fig. 9(b), where upon approaching the cut-off frequency of each flexural wave bandgap, the wave number computed using the homogenization method converges towards the imaginary zero, while the reference solution converges towards the real zero. This deviation in convergence behaviour essentially results from neglecting the rotational inertia, which can contribute to the 2nd-order term in Eq. (66b). This very small detail is, however, irrelevant for most practical analyses.

Computational costs
The computational costs of the ω(k) and k(ω) dispersion analyses using either the homogenized enriched continuum formulation or the standard Bloch approach are shown in Table 4. Note, that for the homogenizationbased analysis only the on-line costs of solving the dispersion problem for the reduced continuum are shown; the off-line simulation costs are disregarded, since this needs to be done only once for a given unit cell. The results demonstrate a significant computational cost reduction provided by the homogenization method for the subwavelength dispersion analysis, with an on-line computational speed that is two orders faster than the standard Bloch analysis. This essentially results from the significantly fewer DOFs in the homogenized reduced-order unit cell model, which governs both the solution and storage costs.

Numerical simulations on a finite LRAM panel
The dispersion analysis presented above relies on the assumption of infinite panel. However, in any practical application, LRAM panels will have finite dimensions, together with boundary conditions, where a Bloch analysis cannot be applied any more. Hence, further studies on the performance of the homogenization method are demonstrated on a finite LRAM panel in steady-state and transient conditions.
The finite LRAM panel sketched in Fig. 10(a) is considered, constructed by the same unit cell as shown in Fig. 6(b), with the material and geometrical parameters given in Tables 1 and 2, respectively. The effective material Table 4 The on-line computational costs of the dispersion analyses on the infinitely long LRAM panel. "N. Points" indicates the number of wave number and frequency points used in the ω(k) and k(ω) forms, respectively. "Deviation" gives the maximum deviation of the homogenization method in the phase velocity at the highest considered frequency 2000 Hz.  properties are thus the same as those presented in Table 3. Only modes n = 2, 3, 5, 6 will be included in the subsequent analysis, since the other (rotational) modes do not couple with the flexural and bulk waves, as can be concluded from the almost zero effective coefficients of modes 1 and 4 in Table 3. The finite LRAM panel consists of 50 unit cells, which results in a panel length L tot = 500 mm. The panel has one unit cell in the thickness direction, with a thickness H = 8 mm. The length to thickness ratio of the panel is therefore 62.5, which is considered to be sufficiently large to justify the thin beam/shell assumptions. The top and bottom surfaces of the panel are free. At the right end of the panel, a hinge constraint is applied in the middle, while at the left edge the displacement in the transverse⃗ e c 3 -direction is prescribed, see Fig. 10(a). In the⃗ e c 2 -direction, the panel is assumed to be infinitely long; the problem can thus be reduced to a beam at the macroscale, while for the microscopic unit cell model, the plane strain condition has to be assumed in that direction.
The unit cell finite element mesh is the same as shown in Fig. 6(b). For the DNS, which will be used as the reference solution, this results in a total number of 71188 DOFs, using a standard 2D plane strain FEM formulation.
For the macroscale beam problem, the weak form of the momentum balance is obtained from Eqs. (11) and (12), with account for the kinematic relation (4), and disregarding the rotational inertia: ∫ which should hold for all admissible, C 1 -continuous test functions δ ⃗ u M,r . The primary unknown kinematic field in Eq. (67) is the displacement field ⃗ u M,r , which should also be at least C 1 -continuous. For the spatial discretization of the macroscopic beam, the isogeometric (IGA) approach [56] has been employed here. The higher-order continuity properties of IGA formulations are especially advantageous for modelling structural elements, e.g. beams and shells [57], and vibration and wave propagation analyses [58]. The macroscale beam problem under consideration has been discretized using 16 1D spline elements with 3rd-order B-spline basis functions. For convenience of the boundary condition implementation, the control net is chosen in a specific way such that the displacement constraints can be directly inserted in the displacement coefficient column. A convergence analysis (not presented here) with respect to the number of spline elements has been carried out based on the eigenfrequency analysis of the beam with free boundaries.
The internal dynamics evolution Eq. (58) for the enrichment variables ζ n are ODEs coupled to the macroscopic momentum balance PDE. This coupled PDE-ODE system of equations can be solved in several different ways, e.g. by treating the enrichment variables ζ n as either the internal variables evaluated at the integration points, or as spatially discretized primary unknown fields (see [59] for the detailed comparison of these two methods in the context of an enriched diffusion formulation). Here, the latter approach has been adopted, whereby the enriched field variables ζ n have been discretized using the same isogeometric 3rd-order B-spline basis functions (although this high order of continuity is not required for ζ n , this choice has been made here for convenience of implementation).

Steady-state frequency response analysis
First, the steady-state response of the finite LRAM panel under harmonic excitation is analysed. In this case, the displacement and force fields are harmonic in time, and thus the problem can be transformed to the frequency domain. The considered frequency range is 0-2500 Hz. The locations at 2/5L tot and 4/5L tot (i.e. at the locations of the 20th and 40th unit cells) are selected as the "gauge" points for the response analysis.
The macroscopic frequency response functions computed from the homogenized enriched beam model at these two points are shown in Fig. 11 and compared against the full-scale reference DNS solutions. To enable comparison, the DNS results have been taken as the corresponding through-thickness averages in the matrix material. Many global structural resonances and anti-resonances are observed in Fig. 11, revealing a high complexity in the frequency response characteristics of the finite structure for practical LRAM panel design. As expected, two wide stop bands emerge, that exhibit a strong attenuation effect on the flexural wave. Their frequency positions and widths adequately match those predicted in the dispersion analysis, see Fig. 8 or 9 .
Furthermore, it can be observed that the homogenization method provides highly accurate results that are close to the reference DNS solutions, both, inside, as well as outside, the frequency bandgaps. The discrepancy between the two solutions shows a minor increase with frequency, which can be attributed to the small deviation in phase velocity, as discussed in Section 7.2.3. close to the applicability limit of the homogenization model for the considered LRAM unit cell based on the separation of scales requirement and the number of modes included in the dynamic basis (all below 2000 Hz). Indeed, beyond 2000 Hz the deviation increases faster and reaches approximately 3.2% at 2500 Hz.
The steady-state macroscopic transverse displacement distributions along the LRAM panel at six representative excitation frequencies 250 Hz, 500 Hz, 1000 Hz, 1750 Hz, 2000 Hz and 2500 Hz obtained through the homogenization method and DNS are shown in Fig. 12(a). For visual reference, the deformed contours of all unit cells constituting the LRAM panel from the DNS are also shown in Fig. 12(b). At the frequencies 250 Hz, 1000 Hz, 2000 Hz and 2500 Hz the local resonance effects are not active and the LRAM panel essentially behaves like a host panel. The structural deformation is clearly dominated by the flexural waves, with their wavelengths rapidly decreasing with increasing frequency. The panel undergoes significant structural distortions, with actual magnitude comparable to the unit cell dimensions, especially in the vicinity of the global resonance frequencies.
Another notable characteristic is that in the absence of local resonance, the steady-state macroscopic transverse displacement distribution along the LRAM panel has an almost non-decaying harmonic form. At the excitation frequency 500 Hz, which belongs to the first stop band, a strong spatial decay can be observed, where the flexural wave is significantly attenuated over a very limited distance, while the remaining unit cells are nearly at rest. Such a fast decay is also observed at the excitation frequency 1750 Hz, belonging to the second stop band. Finally, the accuracy of the homogenized model is again emphasized. As discussed before, the (small) discrepancy in the predicted displacements increases with frequency, and is more pronounced at the higher-frequency global structural resonances, aggravated by the frequency shift (see Fig. 11).

Transient analysis
Next, the performance of the homogenization method for solving problems involving general, not purely harmonic, excitations, is demonstrated. In this case, a frequency domain formulation cannot be used in general and the set of governing equations has to be integrated in time. Indeed, for the developed homogenization method the time domain solution procedure can be easily established, since the homogenized effective material properties are not frequency dependent in the present approach. For the time integration, the Newmark method [60] has been used here, with typical parameters α = 0.2, β = 0.36 and γ = 0.7.
The same finite LRAM panel, with the same material properties and boundary conditions as before (see Fig. 10), is considered. The applied transverse displacement at the left edge of the panel u in (t) is given by where U in is the constant amplitude, f is the frequency and T is the period. This function is shown in Fig. 13. The applied displacement excitation, therefore, consists of a constant term and a harmonic component and is applied during a limited duration of four periods. The transient analyses using the homogenized enriched beam model and the fully resolved panel model have been performed at several excitation frequencies f = 250 Hz, 500 Hz, 1000 Hz, 1750 Hz, 2000 Hz and 2500 Hz. The resulting transverse displacement time histories at two gauge points, 2/5L tot and 4/5L tot , respectively, are shown in Fig. 14. The displacement distributions along the LRAM panel at t = 4T are presented in Fig. 15(a), accompanied by the deformed contours of all unit cells constituting the LRAM panel from the DNS in Fig. 15(b).
In Fig. 14, the transient nature of the considered problem is clearly visible. The presence of the constant component in the displacement excitation (see Eq. (68)) results in transient response profiles that are quite different from the harmonic forms. In spite of this, the frequencies in the stopping bands (500 Hz and 1750 Hz) are clearly attenuated, as expected due to the local resonance effects. These results are confirmed by Fig. 15, showing significant panel distortions for the frequencies outside the frequency bandgaps and the spatial decay of the wave within the frequency bandgaps of the considered LRAM panel.
Finally, it is emphasized that also in the transient case a remarkably good match is observed between the homogenization-based results and the reference DNS results, see Figs. 14 and 15. Similar to the steady-state case, the match between both methods is excellent at the considered (low) frequencies, while the deviation increases slightly with increasing frequency, which is due to the minor deviation on the phase velocity predicted by the homogenization   method. The deviation is, however, limited. It stays below 1.62% for the frequencies up to 2000 Hz and increases to 3.22% at 2500 Hz, which is, in fact, outside the range of applicability of the homogenization method for the considered LRAM panel.

Computational costs
An indication is given below of the computational cost reduction provided by the homogenization method compared to the full-scale DNS. For the above steady-state and transient simulations on a finite LRAM panel, the IGA discretized homogenized model has been implemented in MATLAB, while for the full-scale FEM-based DNS, COMSOL Multiphysics was used. Therefore, a quantitative comparison of the computational costs cannot be done and the numbers are indicative only. All simulations have been performed on the same platform and no parallel computation was used. The representative numbers resulting from the averaging among six tests of the same type are given in Table 5. As before, for the homogenized method only the on-line time is shown. The numbers demonstrate the significant benefit of using the homogenization method for macroscopic analyses of a finite LRAM panel, with an on-line computational gain up to three orders of magnitude compared to the full-scale analysis. The  homogenization method is expected to be even more computationally attractive for the analysis of larger LRAM panels, especially in the transient regime, when also memory allocation and storage costs will increase significantly.

Conclusions
This contribution presented a computational homogenization method enabling an efficient and accurate analysis of wave propagation problems in LRAM panels. The approach is based on a relaxed separation of scales principle, which restricts the methodology to phenomena governed by localized resonance. The macroscopic LRAM panel is modelled as a thin continuum beam/shell described by proper structural kinematics and momentum balance relations. At the microscale, a LRAM unit cell is considered with lamina-wise in-plane boundary conditions and zero out-ofplane tractions representing the free top and bottom surfaces of the macroscopic panel. Under the assumption of linear elasticity, in the relaxed separation of scales regime, the microscale unit cell problem can be represented by the superposition of a static and a dynamic problem, enabling model reduction. As a result, the macroscale effective material properties can be computed once and for all for a given unit cell. In addition, a new macroscale evolution equation emerges describing the effect of the microscopic internal dynamics on the macroscopic fields through the introduction of new macroscale enrichment variables, with a clear physical meaning, i.e. the modal amplitudes of localized microscale modes.
The performance of the proposed homogenization methodology is examined on an exemplary LRAM panel, of which the unit cell reflects a classical negative effective mass design. Several types of wave propagation analyses have been performed, namely computation of the dispersion spectra using both the ω(k) and k(ω) dispersion analyses, a steady-state frequency response analysis of a finite panel under harmonic excitation and a fully transient analysis of a finite panel under non-harmonic loading. For all considered cases, within the defined applicability regime, the homogenization method has been shown to provide highly accurate results compared to the reference DNS solutions, at a fraction of the latter's computational cost. The presented case studies also display a significant potential of using LRAM panel design for attenuating low-frequency flexural waves, which can cause considerable distortions in thin structures, resulting in undesired acoustic coupling effects or even undermining the structural integrity.
The current numerical examples were restricted to a LRAM beam-like panel, with only one direction of wave propagation. In the future, more general shell-like examples will be considered. Furthermore, more complex metamaterial unit cell designs will be investigated, e.g. non-symmetric unit cell configurations with respect to the mid-plane and unit cells exhibiting negative effective stiffness or double-negative effects. Besides, the current framework can be extended to thick LRAM panels, where the transverse shear effects are non-negligible, accompanied by a more pronounced couplings between the flexural (and torsional) waves and the local resonance.
To summarize, the developed homogenization method leads to a reduced, computationally efficient and accurate, equivalent macroscale model of a locally resonant beam/shell medium in the subwavelength range, thereby providing a considerable advantage for fast analyses and design of LRAM panels for practical applications.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Here it has been taken into account that ∆⃗ X + m = ∆⃗ X − m for each lamina, andĤ M =F M −Î, denoting the in-plane displacement gradient tensor.
In the following, the above constraints will be specified for the case of a 2D plane strain LRAM unit cell, which was used for the homogenization towards a macroscopic continuum beam in the numerical examples presented in this paper. A more general case of a 3D unit cell, representative of a macroscopic plate or shell can also be similarly derived.
A 2D unit cell representing a macroscopic beam is sketched in Fig. A.16. Note, that the internal resonator structure is not shown in the figure for clarity.
From Fig. A.16, it can be readily concluded that for the considered unit cell, the following geometrical relations hold between the in-plane position vectors of the points on the left and right boundaries (marked by the superscripts "L" and "R", respectively) Next, conditions (A.4) are reformulated in terms of the displacements of the "prescribed" nodes ⃗ u m,p . To this end, the corner nodes 1, 2, 3 and 4 (see Fig. A.16) are selected to be prescribed. The in-plane displacements of these nodes are thus prescribed according to Eq. (38). The out-of-plane displacement of the corner node 1 is fixed, while for the nodes 2, 3 and 4 the out-of-plane displacements are left free in accordance with the zero traction boundary conditions (27)  • O 1 : These terms are identically zero, sinceÎ m,p and η m,pÎ m,p belong to the null space set ofB uu m . The former implies the rigid-body linear motion, while the latter one corresponds to the rigid-body rotation.
• O 2 : These terms are generally non-zero but vanish when the unit cell exhibits both in-plane and out-of-plane symmetries, with ⃗ X m,r being its inertial centre. In this case, the reference surface also behaves like a neutral stress surface. The symmetric unit cell design will be considered for the numerical examples presented in this paper and thus these terms will be omitted. Using definitions of h N n , h M n , ⃗ j P n and ⃗ j Q n given in Eqs. (53g)-(53j), allows rewriting (B.9) in the compact form (58).