Frequency domain boundary value problem analyses of acoustic metamaterials described by an emergent generalized continuum

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers. Link to publication


Introduction
Acoustic metamaterials are specially designed composites that exhibit extraordinary elastic wave manipulation properties going beyond conventional materials, leading to many potential breakthroughs in noise isolation, acoustic cloaking, precision sensing, wave-guiding, etc [1][2][3]. The exotic properties of acoustic metamaterials are a result of two main phenomena, namely Bragg scattering and local resonance [1]. Bragg scattering occurs in periodic composites when the wavelength of the propagating wave is comparable to the lattice constant. Opposed to that, materials with embedded micro-resonators, called local resonance acoustic metamaterials (LRAM), exhibit a strong coupling to the propagating wave in the subwavelength regime, i.e. when the wavelength is much larger than the unit cell size, while the affected frequency is typically around the natural frequency of the resonators.
Due to the large separation of scales present in LRAMs, several papers have proposed modified effective elastody-B V. G. Kouznetsova V.G.Kouznetsova@tue.nl 1 Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, The Netherlands namic theories which capture the effect in the form of dynamic effective (frequency dependent) constitutive models [4][5][6][7][8] or, equivalently, an enriched micro-inertial continuum [9,10]. Formal homogenization theories bridging the effective models to the classical elastodynamics have also been proposed by various authors [11][12][13][14][15][16]. Furthermore, exploiting numerical solution methods (e.g. FEM), a computational homogenization (FE 2 ) framework suitable for simulating LRAM problems involving arbitrary, complex micro-structure geometries, transient macro-scale loading on finite size domains, nonlinear material response etc., have also been proposed [15,17]. Alternatively, restricted to linear elasticity, Sridhar et al. [16] combined the framework of Pham et al. [15] with dynamic model reduction techniques resulting in an emergent enriched micro-inertial homogenized continuum. This framework served as a basis for the semi-analytical analysis of LRAM designs [18] and allowed to construct computationally efficient multiscale techniques for large scale LRAM boundary value problems presenting a viable alternative to direct numerical simulations (DNS) [19].
While most efforts in developing new homogenization frameworks have focused on local resonance phenomena, a homogenized description of the Bragg scattering phenomena is much more challenging and has been attempted in only a few papers. The reduced scale separation can be taken into account in asymptotic homogenization by adding highorder correction terms [13,[20][21][22][23]. However, this approach is still limited to the first (acoustic) dispersion branch and in some cases, the second (first optical) branch. In higher frequency regimes, the classical scale separation breaks down rendering standard approaches inapplicable. However, for a periodic micro-structure, spectral decomposition of the fields into slow and fast wavelengths is still possible by means of the Floquet-Bloch theory [24,25]. Direct application of the theory allows to reduce the problem on an infinite domain to a problem on a single periodic unit cell, leading to computationally efficient approaches for dispersion [26][27][28] and transmission analyses [29][30][31].
The Floquet-Bloch theory therefore provided a basis for formulating general homogenization frameworks applicable in the high frequency Bragg scattering regime. One of the pioneering steps in this direction was made by Willis [32][33][34][35] and subsequently by other authors [36][37][38][39][40]. However, it has been demonstrated in [36,41] that this theoretical formulation still fails in capturing dispersion spectra at high frequencies. The bottleneck was identified in the classical homogenization operator, which is usually defined as the spatial average over the unit cell volume. In the work of Boutin et al. [14] and Craster et al. [42] on high frequency homogenization, it was shown that the macroscopic fields can be defined as large scale modulations of the Floquet-Bloch eigenmode at a given frequency. The homogenized equations emerging from such an approach capture the local elastodynamics in a narrow range around this eigenfrequency, even at very high frequencies. This motivated Sridhar et al. [41] to extend this idea and define a generalized homogenization operator as a weighted volume average with respect to a family of orthogonal projection functions composed of a set of Floquet-Bloch eigenmodes computed at selected Brillouin zone points in the frequency range of interest. This concept, combined with the Floquet-Bloch formalism, gave rise to a generalized homogenization framework for acoustic metamaterials [37,41]. The resulting macro-scale balance equation describes an emergent generalized elastodynamic continuum of the micromorphic type, as defined by Eringen [43]. The framework demonstrated the ability to accurately predict dispersion phenomena in a wide range of frequencies, including Bragg scattering, emergent elastodynamic phenomena like rotational waves and the hybridization of Bragg scattering and local resonance effects [41].
In this paper, this general elastodynamics homogenization framework is exploited towards the numerical solution of macro-scale boundary value problems. This is a novel development since only a few works have investigated homogenized boundary value analyses in the high frequency regime [39,40,44]. For the sake of simplicity and compactness, the present analysis is restricted to the frequency-domain analysis in a 1D domain. The micromorphic weak form is obtained via a Galerkin approximation. For spatial discretization, an isogeometric formulation with (open uniform, non-rational) B-splines is adopted. This provides high-order continuity properties leading to superior numerical dispersion properties over the traditional finite element shape functions [45][46][47][48]. Since the general homogenization framework introduces several generalized kinematic macro-scale fields, translation of the physical boundary condition of the full-scale problem to the boundary conditions for the generalized continuum is not straightforward and also non-unique. A unique boundary condition leading to the most accurate result is one that excites the longest wave eigenmodes of the system at a given frequency. This is determined by means of an eigenvalue analysis of a special high-order impedance matrix in the space of possible solutions, yielding the appropriate macro-scale boundary mode of the micro-morphic fields.
The structure of the paper is as follows. The theoretical aspects, including the setup of the 1D full scale problem, outline of the general homogenization framework and the incorporation of the macro-scale boundary conditions, are presented in Sect. 2. The numerical dispersion analysis of the B-Spline based discretization and the dispersion spectra predicted by the discretized generalized continuum macroscale model are demonstrated in Sect. 3. The numerical frequency-domain boundary value analysis of the homogenized generalized macro-scale continuum and its validation against direct numerical simulations (DNS) are presented in Sect. 4 on an example problem. Finally, the conclusions are given in Sect. 5.
The following notations are used throughout the paper to represent different quantities and operations. The partial derivative with respect to a variable y is denoted as (•) ,y and the nth order derivative is denoted as (•) ,y n . Matrices of any type of quantity are in general denoted by (•) except for a column matrix (vector), which is denoted by (• ). A left superscript is used to specify a sub-matrix index of a matrix.

Homogenization of the 1D elastodynamic problem
This section presents the main theoretical aspects of the homogenization framework and the solution strategy, particularized to the 1D case. First, the full-scale elastodynamic equation and the associated boundary conditions are defined for a periodic linear elastic medium. The general homogenization framework introduced in [41] is then outlined and applied to recover the homogenized generalized macroscale continuum problem. The addition of artificial highorder damping to suppress spurious solutions is discussed next. Finally, the weak form of the macro-scale balance is derived.

Problem definition
The reference 1D full-scale heterogeneous problem is illustrated in Fig. 1. Consider an unbounded spatial domain R, with coordinate x. The balance of momentum for an elastic heterogeneous continuum is given by, where σ and p represent the scalar Cauchy stress and momentum density, respectively, and t the time. The linear elastic constitutive equations are given as follows where u is the displacement; C m (x) and ρ m (x) are the elasticity modulus and mass density, respectively. Furthermore, a periodic material distribution is assumed, which is expressed as where is the characteristic lattice constant. A unit cell is obtained by considering an arbitrary continuous subdomain of length . A boundary value problem is constructed on a finite subdomain Ω of R, consisting of N cells unit cells (fractional unit cells are not considered here) with boundaries at x 0 and x L = x 0 + N cells as shown in Fig. 1. Furthermore, the macro-structure is constructed using symmetric unit cells and is therefore symmetric about its center axis (as shown in the figure). This avoids structural defect modes [49] that may not be captured in the present framework. This is not a severe limitation since acoustic metamaterials are manufactured composites and it is expected that they will be designed to prevent undesired resonance modes.
Let ∂Ω = {x 0 , x L } represent the set of boundary coordinates. The general boundary conditions on ∂Ω can be defined as follows.

Traction boundary condition
where ∂Ω p and ∂Ω s are the set of boundary coordinates with prescribed displacements and prescribed tractions, respectively, with ∂Ω s ∪ ∂Ω p = ∂Ω and ∂Ω s ∩ ∂Ω p = {} and, σ (x s , t) the applied traction at (x s , t) and u (x p , t) the prescribed displacement at (x p , t). The initial conditions are not explicitly stated here since they are not relevant for the regime of interest in this paper.

The general elastodynamic homogenization framework
The general elastodynamic homogenization framework introduced in [41] is outlined in the following, particularized to the present 1D setting. The framework utilizes a spectral definition of scales (illustrated in Fig. 2) based on the characteristic length scale of the problem. For a periodic structure, this is readily identified as the lattice constant . Accordingly, the macro-scale (or slow-scale) variation of a field is characterized by all modulations with wavenumbers k (spatial frequency) lying within the (first) Brillouin zone, i.e. k M = {k | |k| ≤ π }, while the micro-scale (or fast-scale) variation is characterized by fluctuations with complimentary wavenumbers k m = {k | |k| > π/ k = 0}. More specifically, the micro-scale variations are spanned by Fourier functions defined on the unit cell domain, and hence parameterized by discrete wavenumbers, k m = n 2π , ∀ n ∈ Z (where Z is the space of integers). Since the two wavenumber sets are exclusive (with the exception of the zero wavenumber), a distinct separation of the two scales is obtained. Note, that unlike classical homogenization theories [50], the requirement of large scale separation | k M |<< π/ on the macro-scale variations is not imposed here.
The above introduced concept of spectral scales can be formalized by the Floquet-Bloch transform [36]. Given a function g(x) ∈ L 2 (R), where L 2 (R) is the space of square integrable functions parameterized on R, the Floquet-Bloch transform F B , of g(x) is given bŷ where notation (•) is used to denote the Floquet-Bloch transform of a function; i denotes the imaginary unit. The function g(x, k M ) satisfies the following periodicity conditions, The corresponding inverse Floquet-Bloch transform (synthesis) is given as follows, Following the periodicity property (8a),ĝ(x, k M ) per definition characterizes the micro-scale variations of g(x). Similarly, the harmonic modulation e ik M x in Eq. (9) describes the macro-scale variations of g(x), see Fig. 2 for an illustration. Hence, the Floquet-Bloch expansion (9) provides a natural spectral decomposition of the micro-and macro-scale variations of a function. Using this spectral scale separation, macro-and microscale functions can formally be defined as follows. A macroscale function, F M (x) is one whose Floquet-Bloch transform is constant with respect to x, i.e.
where δ(k M ) represents the Dirac delta function. Alternatively, a micro-scale function can also be defined as one that is periodic with respect to . Note that the material parameters C m (x) and ρ m (x) introduced previously are micro-scale functions depending on the local coordinate.
Homogenization can now be defined as an operation that returns a macro-scale function. To this end, the Floquet-Bloch average, or FB-average, is introduced based on (9). The FB-average of a function g(x) is given as follows, Hereĝ(x, k M ) is averaged over the unit cell domain resulting in a spatially constant function, which upon taking the inverse Floquet-Bloch transform yields a macro-scale function G M . Different realizations of the macro-scale function can be achieved by taking a weighted FB-average in Eq. (12) with respect to some micro-scale weighting/projection func- This is an important aspect from a physical point of view since by formulating the projection functions based on the characteristic micro-dynamic modes of the system (i.e. the Floquet-Bloch eigenmodes) various emergent elastodynamic phenomena occurring in different frequency regimes can be effectively identified at the macro-scale. This concept is exploited here towards a general homogenization operator as will be defined in the following. At this stage, two important consequences of the definition of a micro-and macro-scale functions and the FB average needs to be emphasized. First, from Eqs. (11) and (13), one obtains that the weighted FB- (14) Second, the weighted FB average of a macro-scale function returns the same function times a constant (depending on the selected weighting functions): from Eqs. (10), (13) and (14), the FB-average of the product of an arbitrary function g(x) and a macro-scale function F M (x) can be expressed as follows, Let φ m be a set of N dof orthonormalized micro-scale functions parameterized by a set P = {1, 2, . . . N dof }. The orthonormality condition can be expressed as where I is the scalar identity matrix of size N dof and (•) T denotes the transpose of a (column) matrix. A general homogenization operator with respect to φ m , S P , is defined based on the FB average (13) as follows, Therefore, the application of S P on a full-scale field f (x) yields a column with N dof generalized macro-scale fields F M . The projection functions φ m span the subspace of all possible micro-scale ( -periodic) functions. A larger P increases the projection space of the homogenization operator allowing more micro-scale information to be captured at the macro-scale. In the limit when P → P com (where P com = {1, 2, . . . ∞}) an invertible mapping is obtained that, when applied on F M , will recover all information of the full scale P com is the inverse of S P com , termed localization operator. However, a complete mapping is neither practical nor desirable (nor possible in most cases) as the number of the required projection functions would need to be too large (infinite in the limit). Thus, an appropriate formulation of φ m is critical in capturing the relevant physical processes occurring in the regime of interest. Indeed, for computational efficiency, one should aim for the smallest possible set of projection functions that still captures the desired physical phenomena with sufficient accuracy. Thus, the precise selection of the projection functions matters. This aspect will be detailed in Sect. 3.3.
As explained above, a unique localization operator, S −1 P cannot be realized since P is incomplete in general. Thus a generalized localization operator, denoted as S −1 Q will be introduced such that, Here Q is a set (not necessarily equal to P) that parameterizes the micro-scale functions used to construct S −1 Q . The localization functions must be constructed such that it adequately captures all the micro-mechanical/dynamic response of the microstructure in the frequency range of interest. In practice, an optimal formulation of φ m as basis functions for homogenization is usually sub-optimal for constructing S −1 Q . A more sophisticated approach for formulating S −1 Q is therefore needed. In [41], the localization operation was expressed as a high-order gradient spatio-temporal expansion with respect to the macro-scale quantities where the individual microscale functions were determined using an energy equivalence between the scales. Here, the considered 1D problem is sufficiently simple, allowing to construct the space spanned by φ m from a set of characteristic dynamic eigenmodes of the microstructure. This is sufficiently accurate and recovers the micro-dynamic response in the frequency range of interest. Accordingly, the following localization operator for Q = P is adopted, The above definition is clearly consistent with the identity (18), as can be verified by making use of the orthonormality condition (16).

The homogenized macro-scale boundary value problem
The homogenization framework presented in the previous section is now applied to obtain the homogenized macro-scale equations for the 1D problem outlined in Sect. 2.1. Since the Floquet-Bloch transform is defined on an unbounded domain, the homogenized equations will be formulated on R first. The application of the boundary conditions will be discussed later. The homogenized balance of momentum is obtained by applying the homogenization operator (17) to Eq. (1), which by using the Gauss divergence theorem yields (see [41] for more details), where Here Σ M , χ M and q M are the columns (of size N dof ) of the generalized macro-scale stress, internal stress and momentum fields, respectively. Next, the homogenized constitutive equations are derived through the sequential application of the localization and homogenization operators. First, the full scale displacement is approximated using the localization operator (19) where * u(x, t) represents the approximation of u(x, t) and U M (x, t) represents the column of N dof generalized macroscale degrees of freedom defined via Eq. (17) as Substituting Eq. (22) into Eq. (2) and applying Eqs. (21) [while observing rule (15)] yields, represent the macro-scale constitutive matrices of size N dof × N dof . Note that the constitutive law is local in nature and the parameters are constants with respect to space and time. Equations (20) and (24) constitute the set of generalized macro-scale continuum equations. Substituting the constitutive Eqs. (24) into (20) yields the governing partial differential equation in terms of the generalized displace- The finite domain Ω containing the homogenized generalized continuum, is considered next for setting up the boundary value problem. The full-scale problem boundary conditions (5) and (6) have to be converted to expressions in terms of the homogenized generalized fields. To this end, the localization expressions of the full-scale displacement and stress are needed. The localized full-scale displacement field was given in Eq. (22). Next, utilizing Eq. (19) the localization of the full scale stress is obtained as * where * σ (x, t) represents the approximation of σ (x, t). Alternatively, the localization of the full scale stress can also be obtained by substituting Eq. (22) The choice between the two approximations, Eqs. (28) or (29) is arbitrary. Equation (28) is here adopted as the approximation for the full-scale stress. With (22) and (28), the full scale boundary conditions (5) and (6) can be expressed in terms of the homogenized fields as follows Macro-scale traction boundary condition Macro-scale displacement boundary condition providing 2 boundary conditions. However, the homogenized problem of the generalized continuum requires in total 2N dof conditions which (assuming N dof > 1) cannot be uniquely determined from the full-scale problem without additional information. Thus, additional conditions must be formulated in order to overcome this intrinsic non-uniqueness. An additional constraint could be the requirement for power consistency at the boundaries. This was used in [40], replacing the traction boundary condition that could not be satisfied due to the occurring boundary layer effects. Alternatively, an additional condition based on the nature of excited solutions can be formulated for determining macro-scale boundary values. This approach will be elaborated in the subsequent section. The boundary power consistency will be analyzed a posterori (see Sect. 4).

Plane wave expansion of the homogenized equations: towards an approach for spectral decomposition and solution filtering
In addition to the physical solution, the homogenized generalized macro-scale continuum also admits several duplicate/spurious solutions. This is due to the definition of the homogenization operation, which introduces multiple macro-scale fields representing a single fine-scale field, resulting in an increased number of admissible wave solutions of the problem. The spurious solutions tend to introduce a numerical error in the solution of the boundary value problem and must therefore be eliminated or attenuated. The decomposition of the frequency-domain solution into the primary and spurious solutions can be obtained through the analytical plane wave expansion of the homogenized equations. This is used to shed light on the nature of the spurious solutions and to formulate an approach to mitigate these solutions in a numerical setting.
The plane wave expansion of U M at frequency ω, can be expressed as where rV M (ω), r k e M (ω) and r y M (ω) are the eigenmode (polarization), eigen-wavenumber and amplitude of the r th wave (r = 1, . . . , 2N dof ), respectively. The corresponding eigenvalue problem is obtained by substituting Eq. (32) into Eq. (26) yielding, Here ω is treated as a parameter when solving the above eigenvalue problem. The amplitudes r y M are to be determined such that the solution (32) satisfies the boundary conditions (30) and (31). However, as discussed in the previous section, the solution is not unique for N dof > 1 since the number of unknowns is greater than the number of boundary conditions given by the full-scale system. This implies that not all wave modes in the expansion are relevant for the solution. In fact, it can be shown that 2 solutions exist whose wavenumbers (at least the real part) lie within the Brillouin zone, while the rest lie beyond it [27]. Thus, based on the definition of a macro-scale function as presented in Sect. 2, it is clear that the solutions outside the Brillouin zone constitute the spurious contribution and must therefore be discarded. However, a priori identification of these solutions via plane wave expansion is only possible for simple boundary value problems. A more general approach is therefore needed in order to adapt the method towards complex boundary value problems involving 2D/3D macro-structures. It is therefore desirable to modify the system of equations such that it retains the primary wave solutions while the spurious solutions are all made evanescent. To this end an artificial forcing term proportional to the high-order spatial derivatives of the macro-scale velocities is added to Eq. (33) as follows where the first term has been added with ζ , γ and Z M as parameters. For γ ≥ 2, the above term provides a highorder damping that is proportional to the (2γ )th power of the wavenumber, thus high wavenumber solutions will exhibit a large evanescent component, thus being effectively attenuated within the bulk of the system. The higher the value of γ , the stronger the bias towards higher wavenumbers. The choice of the coefficient Z M is arbitrary but a suitable formulation that provides proportional damping on all generalized macro-scale displacements was found to be the homogenized acoustic impedance matrix The square root in the above equation represents a matrix square root (i.e. B = A ⇒ A = B B ). The parameter ζ is tuned for a given numerical discretization such that the primary solution remains unchanged, while the spurious solutions are sufficiently attenuated. (Explored further in Sect. 3).
Introducing the higher-order damping into the macroscale balance Eq. (20) leads to where the subscript (•) ,x 2γ denotes the 2γ th derivative w.r.t x. The addition of the high-order damping term increases the order of the macro-scale equation from 2 to 2γ . Since the higher-order term has only been introduced to provide artificial damping, a natural assumption on the higher-order boundary conditions is to require all higher order spatial derivatives of the macro-scale generalized displacements at the boundary to vanish, i.e.

Weak form of the generalized macro-scale balance
Finally, the weak form of the homogenized macro-scale equations facilitating the numerical discretization of the problem is obtained as follows. First, integrating Eq. (36) over Ω with respect to the test function δU M gives, The weak form is then obtained by applying the Gauss divergence theorem and taking into account condition (37) giving, For the purpose of the subsequent dispersion analysis, the weak form on the unbounded domain R is obtained by replacing Ω with R and dropping the boundary terms (since all the macro-scale fields are assumed to be zero at ±∞) in the above expression. Accordingly one obtains,

Numerical dispersion analysis
This section discusses the numerical dispersion analysis of the macro-scale equations. The purpose of this analysis is two-fold, namely to study the dispersion characteristics of the homogenized discrete problem versus the homogenized continuum formulation and versus the reference Bloch analysis on the unit cell. First, the discrete balance is obtained on an unbounded domain. The numerical dispersion equation is then derived by applying the discrete Fourier transform on the discretized balance equation. Next, the analytical spectral characteristics of various order B-spline approximations are investigated. Finally, the numerical dispersion spectrum is validated against the direct Bloch analysis for an example unit cell. The influence of the high-order damping term is also illustrated.

Discrete homogenized continuum model on an infinite domain
The generalized displacement field U M on R is discretized as follows, Here p (x) represents the pth order B-spline polynomial shape function with a compact span Ω span = {x ∈ R | −( p + 1)h/2 ≤ x ≤ ( p + 1)h/2}, h is the spacing between the control points and (n) U d M are the discrete amplitudes associated to the nth shape function. The B-spline shape functions (for p = 2) on an unbounded domain are illustrated in Fig. 3. The explicit expressions of the B-spline functions can be found in literature, e.g. [48,51].
Substituting Eq. (41) into the weak form (40) in combination with the homogenized constitutive Eq. (24) and using a Galerkin approximation for the test functions yields the following discretized stencil equations, Note that the summation over the integer space in Eq. (41) reduces to a summation over the span of one shape function in Eq. (42). Furthermore, the evaluation of (n) D requires shape functions with a minimum polynomial order p = γ . Since γ > 2, linear shape functions cannot be used to capture the high-order damping. The integration is performed numerically using a standard Gauss integration scheme corresponding to the polynomial order of the given shape function, where the elements were defined as the domain between two control points, thus Ω span contains p + 1 elements.

Spectral characteristics of the discretization
The discrete Fourier expansion of (n) U d M is given as follows whereÛ d M are the coefficients of the discrete Fourier transform of (n) U d M . Substituting the above expression into Eq. (42) and gathering the terms with respect to the Fourier coefficient gives, are the functions characterizing the spectral behavior of a given B-spline discretization for the homogenized generalized continuum. Comparing the discrete problem (45) to the equivalent continuum Eq. (34) reveals that the numerical discretization accurately approximates the continuum dispersion spectra in the regime where First, the functions κ 1 and κ 2 are computed for various p refinements and plotted against k M in Fig. 4a, b. Figure 4a, b, reveal that the characteristic functions of the discretized problem recover the continuum spectra in the low wavelength regime k M << π/h and degrade in accuracy for k M approaching the discretization limit (π/h). Furthermore, the overall accuracy significantly improves if polynomial order is increased. However, κ 1 fails to converge at k M = π/h regardless of the refinement order p. This is a direct consequence of discretizing equations with odd-order spatial derivatives. Under Galerkin projection, the discretized displacement field becomes orthogonal to its odd-order approximation in the limit k M → π/h. Thus the numerical dispersion spectrum cannot converge to the actual spectrum at wavenumbers approaching the discretization limit. As discussed in Sect. 2.4, the physical macro-scale solution lies within the Brillouin zone with wavenumbers |k M | ≤ π/ . Therefore, the h refinement should be significantly smaller than in order to avoid errors due to numerical dispersion affecting the solutions within the Brillouin zone. Increasing the order p of the B-spline approximation improves the dispersion characteristics and would therefore allow for a relatively coarse discretization.
Next, κ 2γ is computed for the parameter combinations p = 2 and γ = 2, p = 3 and γ = 2, and p = 3 and γ = 3 and plotted in Fig. 4c. The response of κ 2γ is only sensitive to high wavenumbers. As expected, the sensitivity improves significantly for higher γ (with p kept constant). Similarly, a p refinement also increases the sensitivity of the damping to higher wavenumbers, although the effect is less pronounced.

Validation against Bloch analysis
The numerical dispersion spectrum is now validated against Bloch analysis for an example unit cell model. An ω − k analysis is performed where the dispersion eigenvalue problem is solved by treating k as the eigenvalue and ω as the The dispersion analysis is carried out on the example unit cell micro-structure shown in Fig. 5, i.e. a two-phase composite with a hard heavy phase and a compliant light phase. The choice of the particular unit cell configurations does not play a role in a dispersion analysis but it is given for the sake of completeness (as it will matter in a boundary value problem). The unit cell in configuration A (Fig. 5) was selected here for the analysis. The values of the material and geometrical parameters are described in Table 1.
The reference Bloch analysis is carried out using the FEM approach described in [27]. The weak form of the dispersion equation is given as follows where w represents a test function, (•) * denotes the complex conjugate and T is the unit cell domain. The unit cell is discretized using standard linear finite element shape functions with a uniform element size of 0.01 , which ensures adequate convergence of the computed spectra within the Brillouin zone and in the analyzed frequency range of 0-5000 Hz.
The projection functions φ m are formulated using the Floquet-Bloch eigenmodes at the symmetry points and X (i.e./ k M = 0 and k M = π/ respectively) computed using the same FEM unit cell model as for the reference Bloch analysis.
soft phase hard phase A.

Fig. 5
The considered 1D unit cell represented in two possible symmetric configurations Length (mm) 10 The first 5 modes and the real components of the first 4X modes are shown in Fig. 6. After orthonormalization of these functions using the Gram-Schmidt procedure enforcing condition (16), 9 projection functions result, corresponding to 9 generalized degrees of freedom. The analysis is first carried out on the homogenized model without artificial damping. The ω − k dispersion spectra of (i) the considered unit cell computed using the Bloch analysis, (ii) the homogenized continuum dispersion Eq. (33) and (iii) the homogenized numerical dispersion Eq. (47) using cubic ( p = 3) B-splines with h = 1/3 and ζ = 0 are shown in Fig. 7. The homogenized continuum dispersion spectrum matches exactly at all frequencies, at least with respect to the primary solutions (i.e. within the Brillioun zone). The homogenized numerical spectrum also predicts the primary solution accurately within the Brillioun zone, thereby validating the numerical dispersion model for the given discretization parameters. Additional spuri- ous solutions exhibited by the numerical spectrum at high wavenumbers still emerge as a consequence of the discretization approximation. The influence of the artificial high-order damping on the numerical dispersion model is assessed in Fig. 8. A 6th order damping term is adopted i.e. γ = 3 with ζ = 0.00002. The selective damping of the spurious solutions is clearly visible, while the solutions within the Brillioun zone remain unaffected. For the given parameters, most solutions beyond the Brillioun zone exhibit a spatial evanescense of magnitude greater than 0.05π/ . Thus, the artificial higher order damping will attenuate the (residual) spurious short wave solutions (complimenting the least impedance boundary

Discrete enriched continuum model on a finite domain
The pth order B-spline approximation (illustrated in Fig. 9 for cubic B-splines) of U M defined on Ω is expressed as follows, where f p (x) is a column containing N D = N cells /h + p, pth order B-Spline shape functions defined on Ω, I the identity matrix of size N dof , U d M , is a column containing N Gdof = N dof N D discretized degrees of freedom and ⊗, the Kronecker product operator. Note the matrix notation adopted here in the finite domain case for the sake of convenience as opposed to the stencil notation used in Sect. 3 (Eq. 41). A uniform control point distribution with spacing h is assumed except at the boundaries where the spacings are collapsed to recover the delta property, i.e. 1 f p (x 0 ) = N D f p (x L ) = 1, assuming each shape function is labeled in an ascending order from left to right. Substituting Eq. (24) and Eq. (50) into (39) yields the following discrete balance equation for the homogenized enriched continuum on the finite domain Ω where

Frequency domain analysis and validation
A frequency-domain analysis can be performed by assuming the following solution for the macro-scale displacements and tractions where (•) is used to indicate a frequency-domain solution. Substituting Eq. (53) into Eq. (51) yields the following frequency-domain problem, where

Determination of macro-scale boundary conditions via the least impedance boundary mode approach
This section describes an approach that exploits the highorder impedance matrix D (Eq. 52e) for determining the macro-scale boundary conditions that minimize the appearance of spurious solutions at the boundaries. First, the general solution for the boundary displacements at both ends is computed. The solution is then projected onto D and the eigenmodes corresponding to the smallest two eigenvalues are retained. In principle, these two eigenmodes must be the primary long wave solutions since D is least sensitive to them (as discussed in Sect. 2.4 and 3). Hence, these eigenmodes represent the macro-scale boundary modes of interest at a given frequency. The amplitudes corresponding to the two boundary modes are then uniquely determined by applying Eqs. (30) and (31). The least impedance boundary mode approach is described in more detail in the following. Let the degrees of freedom be partitioned into the boundary and internal dofs represented as b (•) and i (•), respec-tively. Let S be a N Gdof ×2N dof matrix representing the solution to Eq. (54) for different applied generalized boundary displacements i.e.
given as follows, where I represents an identity matrix of size 2N dof . A new reduced high-order impedance matrix D red is thereby obtained by projecting S onto D as follows where (•) * represents the complex conjugate. Let λ j and X j be, respectively, the jth eigenvalue and eigenmode of D red obtained from the following eigenvalue problem, The least impedance boundary mode approach states that the macro-scale boundary values should be determined by a linear combination of the eigenmodes of D red corresponding to the smallest two eigenvalues of the system. Assuming that the eigenvalues are sorted in ascending order according to their magnitudes, this implies that the boundary displacements can be expressed as bŨ d where α 1 and α 2 are unknowns determined by applying the boundary conditions (30) and (31). Thus, substituting Eq. (57) and Eq. (63) into Eq. (54) and applying localization at the boundaries one obtains, where the boundary degrees of freedom have been further partitioned into the nodes on which prescribed displacements and tractions are applied, represented by p (•) and s (•), respectively. Determination of α 1 and α 2 from the solution of the system of Eq. (64) yields the solution to the macro-scale boundary value problem written as follows, where the subscript "imp" is used to indicate that the solution has been derived via the least impedance mode approach.

Numerical example
In the following, several numerical analyses are performed using the homogenized model, which are then compared with DNS. First, the frequency domain response of the macrostructure constructed by sequentially joining the unit cells in configuration A (see Fig. 5) is computed. Next, the procedure is repeated for the unit cell configuration B. The following choices apply to both simulations.
• The macro-structure is constructed using N cells = 10 cells. • A unit displacement is prescribed on the left boundary (i.e.ũ 0 = 1) while a traction free condition is assumed on the right boundary (i.e.σ L = 0). The applied frequency ranges from 0 to 5000 (Hz). The least impedance boundary mode approach is applied. The (reaction) full-scale traction at the left boundaryσ 0 and the full scale displacement at the right boundaryũ L are computed during the analyses and the results are plotted in Figs. 10, 11, 12 and 13. Furthermore, the full-scale and homogenized powers at x 0 , P full and P homo respectively, given as are also plotted for the two simulations in Figs. 11 and 13, respectively. This gives an indication of the resulting mismatch in terms of boundary power. Figures 10 and 12, reveal that the (presented) homogenized results in terms ofσ 0 and u L exhibit an overall excellent match with DNS for both unit cell configurations. The responses for the two configurations (b) Fig. 11 The full-scale versus homogenized power for micro-structural configuration A. Note that iω is dropped from the power term in order to make it comparable toσ 0 A and B (Fig. 5) are not identical especially at higher frequencies (as seen in the traction plot), demonstrating the ability of the general homogenization framework to recover the influence of the material phase distribution, unlike other attempts presented in the literature [39,40]. The power plots shown in Figs. 11 and 13 reveal a discrepancy between the homogenized boundary power of the DNS and homogenized model for frequencies beyond the first bandgap. This inconsistency appears to be an inherent feature of the Floquet-Bloch homogenization, as has been previously pointed out in the literature. For example, in the work of Srivastava et.al. [39,40] an inherent inconsistency was also observed in the boundary tractions and displacements for equal powers. The origin of this apparent power inconsistency is so far unresolved and beyond the scope of the present paper. It also brings to question whether the homogenized boundary power (Eq. 67) can be interpreted as its classical counterpart, i.e. ifP homo =P full holds. The accurate simulation of the transmission response of the example metamaterial structure demonstrates the applicability of the general homogenization framework towards the analysis of boundary value problems in the frequency domain. A computational cost benefit analysis with respect to DNS is not attempted here since the simplicity of the 1D problem considered does not allow to realistically assess the full computational efficiency potential. However, a remarkable computational cost reduction may be expected for more complex problems involving 2D/3D micro and macro-structures. (b) Fig. 13 The full-scale versus homogenized power for micro-structural configuration B. Note that iω is dropped from the power term in order to make it comparable toσ 0

Conclusion
This paper presented a numerical frequency-domain boundary value analysis exploiting a general homogenization framework for the analysis of a 1D heterogeneous problem. The main challenge faced here was the emerging nonuniqueness of the boundary conditions corresponding to the full-scale problem. The general homogenization operator extracts a family of macro-scale fields from a full-scale field. This allows accurately modeling of the dispersion phenomena in a wide range of frequencies. However, the resulting homogenized macro-scale equation requires more boundary conditions than those given by the full scale problem.
An additional requirement on the boundary mode to only contribute to/couple with the primary solution and not additional spurious solutions was proposed in order to uniquely determine macro-scale boundary conditions. To this end, a least impedance boundary mode approach was proposed as an additional step to remedy the problem of the boundary conditions non-uniqueness. The approach was based on the concept of the spectral bias of the high-order spatial gradient operator towards short wave solutions. Based on this, a high-order impedance operator has been proposed, which was used to attenuate the spurious short wave solutions both in the bulk domain (in the form of the high-order damping) and at the boundaries (in form of the minimum impedance mode procedure). A Galerkin B-spline discretization was used to approximate the homogenized macro-scale equations. The methodology was compared with DNS for both numerical dispersion and frequency-domain boundary value analysis on an example periodic 2-phase laminate composite structure. The numerical dispersion analysis revealed that cubic B-splines with a refinement of 1/3 of the unit cell size yield a sufficient accuracy within the Brillouin zone of the spectrum. Two different symmetric unit cell realizations for a given material phase distribution were used in the simulation. The homogenized model showed excellent match with the DNS in capturing the transmission response in the frequency range of the analysis for both cases. A distinct response emerged for the two unit cell configurations, thereby illustrating the ability of the homogenization framework in recovering the influence of the material phase distribution at the boundaries. In spite of the encouraging results obtained, the framework is apparently not (yet) energetically consistent. This needs further considerations when applying the homogenization framework to other settings (e.g. transient analysis and 2D/3D problems). This will be addressed in future work. Overall, the results demonstrate a promising potential of the methodology towards boundary value analyses.