An efficient multiscale method for subwavelength transient analysis of acoustic metamaterials

A reduced-order homogenization framework is proposed, providing a macro-scale-enriched continuum model for locally resonant acoustic metamaterials operating in the subwavelength regime, for both time and frequency domain analyses. The homogenized continuum has a non-standard constitutive model, capturing a metamaterial behaviour such as negative effective bulk modulus, negative effective density and Willis coupling. A suitable reduced space is constructed based on the unit cell response in a steady-state regime and the local resonance regime. A frequency domain numerical example demonstrates the efficiency and suitability of the proposed framework. This article is part of the theme issue ‘Current developments in elastic and acoustic metamaterials science (Part 2)’.


Introduction
Locally resonant acoustic metamaterials (LRAMs) exhibit unconventional material behaviour when their response is observed at a macro-scale.These include negative effective mass density (Liu et al. [1]), negative effective bulk modulus (Fang et al. [2]) and double negativity [3,4].Local resonance phenomena at the micro-scale constitute the underlying mechanism that triggers such anomalous behaviour and can be tuned to take place in the subwavelength regime, providing compact designs relative to the target wavelength [5].Typically, when an effective material property becomes negative, acoustic waves cannot propagate, providing an attractive passive solution for wave mitigation, as extensively discussed in Gao et al. [6] with promising applications for sound insulation and sound absorption, e.g.aeronautics and building acoustics.
A popular emerging type of LRAM is the so-called labyrinthine metamaterials that owe their name to the induced coiling up of the wave path, first proposed by Liang & Li [7].The subwavelength character of labyrinthine metamaterials emerges from the space-coiling property, effectively creating a fluid medium with a low effective sound speed relative to the background fluid.When these structures are arranged sparsely, local resonances take place, providing enhanced sound insulation, first demonstrated by Cheng et al. [8].Later, bio-inspired [9] or fractal [10] patterns have been presented, as well as multiple other designs [11][12][13][14][15][16][17], just to mention a few.Labyrinthine metamaterials are particularly interesting because they can be easily manufactured, e.g. using three-dimensional printing, with high porosity, leading to efficient, lightweight sound insulators.Additionally, their fluid-based resonance mechanism takes place in the absence of deforming resonant structures, which removes potential fatigue issues.
Metamaterial wave dispersion properties are basically characterized by considering an infinite lattice composed of periodically repeated unit cells.Imposing the so-called Bloch-Floquet periodic boundary conditions allows one to understand the wave propagation behaviour by solving a problem for a single unit cell, thus, at a low computational cost, e.g.see the finite element implementation in Elford et al. [18].Alternatively, the effective material properties can be back-calculated from one-dimensional setups [19][20][21].However, the characterization in an infinite domain is rather far from realistic applications where metamaterial devices have three-dimensional boundaries, complex loading conditions and interfaces to other components.Yet, a finite-sized domain problem creates a major concern of computational affordability due to an excessively large number of degrees of freedom required to resolve the full metamaterial structure composed of many unit cells in the time and frequency domain.This identified computational challenge calls for efficient computational methods that enable the full exploration of metamaterial functionalities at the engineering scale, as underlined by Palma et al. [22] in the context of integration of metamaterial technology into aeronautic applications.
Dynamic homogenization is a potential solution to these issues and has now achieved a basic level of maturity.For example, the elastodynamic homogenization theory initiated by Willis [23] and extensively revised by Nassar et al. [24] and Sridhar et al. [25] proves to be exact for an infinite continuum.It defines a macroscopic medium providing macroscopic non-standard constitutive laws for heterogeneous materials.Nonetheless, effective properties in the infinite domain are not necessarily directly transferred to the finite domain, as analysed in Srivastava & Nemat-Nasser [26] and Sridhar et al. [27].In Muhlestein et al. [28] and Sieck et al. [29], Willis homogenization results from elastodynamics are specific for fluid acoustics, where pressure is a more natural representation of the waves.
Several homogenization methods have been proposed for LRAMs composed of Helmholtz resonators in which the fluid is the wave carrier.Hu et al. [30] developed a macroscopic model by deriving analytical formulas of effective properties using a simple cylindrical tree-component geometry unit cell based on a two-step homogenization scheme.Another approach is the rigorous asymptotic homogenization based on scale separation, which has been used to describe wave propagation in rigid porous media with embedded Helmholtz resonators by Boutin [31], and later validated experimentally in Boutin & Becot [32].Based on ensemble averages, Lafarge & Nemati [33] developed a homogenization method based on acoustics-electromagnetics analogy, providing a rigorous macroscopic description of rigid porous media.The method was capable of describing sound propagation in a sequence of Helmholtz resonators by Nemati et al. [34] and in an array of rods in the Bragg scattering regime, i.e. beyond the long wavelength limit [35].These insightful and physically sound homogenization 2 royalsocietypublishing.org/journal/rsta Phil.Trans.R. Soc.A 382: 20230368 approaches mainly focus on determining frequency-dependent and wavenumber-dependent material properties of the corresponding macro-scale constitutive models.However, their application to solving finite-size metamaterial problems has not been presented.
Reduced-order computational homogenization, e.g.Sridhar et al. [36], is a multiscale technique for linear problems that provides an equivalent macro-scale continuum with non-standard constitutive relations retrieved from microstructural unit cell simulations.It belongs to a class of homogenization methods based on the definition of a representative volume element, as extensively discussed in Geers et al. [37,38], and it was first applied to solid metamaterials by Pham et al. [39] employing the FE 2 computational homogenization approach.Reduced-order homogenization allows overcoming computational cost issues, leading to an effective continuum for LRAMs by appropriately transferring the micro-scale inertia to the macro-scale model, which can be straightforwardly implemented, e.g. using the finite element method.Reduced-order computational homogenization provides a macroscopic primary field that approximates the average of the non-resonant domain of the unit cell, and it is based on the so-called relaxed separation of scales, similar to the approach of asymptotic homogenization broadly reviewed in Boutin et al. [40].However, until now, the computational homogenization method has only been applied to elastodynamics, i.e. wave propagation in solid metamaterials described by means of a displacement-based formulation.In contrast, the pressure field naturally describes metamaterials relying on fluid resonance, such as labyrinthine metamaterials, revealing the need for a pressure-based framework.
In this article, we propose an efficient macro-scale model for acoustic metamaterials in which pressure is the primary field, as is typically the case in acoustics.The present work provides an enriched continuum capable of capturing negative effective properties and Willis coupling.It enables assessing the performance of metamaterials beyond the conventional infinite domain restriction for both transient and time-harmonic wave propagation problems.
The article is organized as follows: first, the governing equations of the macro-scale problem are defined in §2.Second, the unit cell governing equations and the model reduction based on subwavelength resonances are elaborated in §3a,b, respectively.Next, the computational homogenization framework is presented by coupling macro-and micro-scales in a variationally consistent manner based on the extended Hill-Mandel condition in §3c, followed by the development of an enriched continuum that describes LRAMs in §3d.Finally, numerical examples show the comparison between the full-field response and the response of the proposed homogenized continuum in infinite and finite domains in §4.Conclusions are given in §5.
For clarity, the mathematical notation is defined in this paragraph.The material derivative is equivalent to the spatial time-derivative due to the use of linearized equations for the pressure fluctuations in the absence of mean flow, see Goldstein [41]; thus, they can be used interchangeably as ∂ ∂t ≡ ˙.The mathematical objects used are scalars (zeroth-order tensor), vectors (first-order tensor) and second-order tensors and are represented as a, a → and A, respectively.
The inner product is denoted as , where δ ij is the Kronecker delta.The inverse and transpose of a second-order tensor are denoted by • −1 and • T , respectively.The notations • ∼ and • are used for column and matrix assem- blies, respectively.The entries of columns and matrices may be scalars, vectors or tensors.The inverse of a matrix is denoted by • −1 .The transpose of a column and a matrix are denoted • ∼ T and • T , respectively.Macroscopic quantities are identified with the subscript 'M', while microscopic quantities do not have such a subscript.

Macro-scale homogenized problem
Consider acoustic wave propagation through a fluid saturating rigid, undeformable channels, a typical problem in acoustics and acoustic metamaterials.The acoustic wave is supported by the fluid, usually air or water.Considering that the numerical simulation of wave propagation through a finite-sized metamaterial, e.g.depicted in figure 1a, has prohibitively high computational costs, the reduced-order homogenization technique developed here aims to propose a computationally efficient equivalent macro-scale problem.The idealized effective problem, figure 1b, is developed based on the analysis of the response of a micro-scale unit cell, figure 1c.The primary field variable for the acoustic metamaterial case is the macro-scale pressure disturbance p M x M → , t (i.e. the pressure disturbances around a quiescent background in the absence of net flow).The relevant macro-scale conservation equation is the acoustic wave equation written in a non-conventional form where ϵ¨M is the second time-derivative of macro-scale volumetric strain, and U → ¨M is the macro- scale acceleration.The constitutive models for the conserved quantity ϵ¨M and its flux U → ¨M will be obtained from the upscaling relations of the micro-scale fields, extracting the essential dynamics.In general forms, macro-scale constitutive relations can be written as functions of the macro-scale pressure disturbance p M x M → , t , pressure gradient ∇ M → p M x M → , t and potential ɴ internal enrichment variables η j x M → , t , with j = 1, . ., ɴ, and their respective time derivatives as The internal enrichment variables must obey ɴ evolution equations (2.3) g j p M , p ˙M, p ¨M, ∇ M → p M , ∇ M → p ˙M, ∇ M → p ¨M, η j , η ˙j, η ¨j = 0, j = 1, . . ., ɴ that are also to be obtained from the micro-scale.The macro-scale initial boundary value problem is described by equations (2.1)-(2.3)subjected to initial and boundary conditions.The boundary conditions are then given in terms of pressure (Dirichlet) or gradient of pressure, i.e. velocity (Neumann) boundary conditions, whereas the initial conditions are given for the pressure field, the enrichment variables and their first-time derivatives.In the case of a homogeneous fluid, the conventional representation of the acoustic wave equation in terms of acoustic pressure is achieved after substituting the relations (2.3) and (2.2) into (2.1).Therefore, the proposed homogenization approach provides an effective fluid behaviour in rigid channels.The core of the proposed framework is to provide a systematic approach to determine the macro-scale constitutive response (2.2), together with the evolution equation (2.3), from the unit cell analysis.It will be achieved through the definition of a scale coupling and homogenization procedure established in §3.

Micro-scale unit cell problem
where ϵ˙ is the volumetric strain rate and U → ¨ the fluid acceleration at the microscopic unit cell scale.These quantities are related to the primary micro-scale unknown field p x → , t , i.e. the fluid pressure disturbance around the background equilibrium (see Landau & Lifshitz [42]), through the constitutive relations of a compressible inviscid fluid where ρ f and K f are the fluid mass density and bulk modulus at the background equilibrium state, respectively.The first equation of (3.2) is the combination of the linearized continuity equation with the linearized state equation, and the second is the linearized balance of linear momentum.Substitution of (3.2) into (3.1)leads to the conventional form of the acoustic wave equation.The dissipation due to viscous effects and heat conduction will be disregarded here.
As pointed out in [43,44], the inviscid response gives the leading-order behaviour of the fluid's effective motion in porous materials for frequencies higher than a minimal frequency , where ν f is the kinematic viscosity of the fluid and d the channel's characteristic width.The sound hard boundary condition is imposed at the interface between the fluid and rigid solid domains as where n → fs is the normal to the interface.Finally, the resulting partial differential equation, (3.1), with (3.2) and (3.3), is complemented by boundary and initial conditions that are to be derived from the scale coupling relations discussed in §3c.

(b) The discretized micro-scale problem projected on a reduced space
For the unit cell, a standard finite element discretization scheme is used for the governing equation (3.1), defining the micro-scale problem with so-far-undefined boundary conditions.The procedure leads to the system of second-order ordinary differential equations written in matrix format as royalsocietypublishing.org/journal/rsta Phil.Trans.R. Soc.A 382: 20230368 where Q is the matrix related to the volumetric strain energy and H the matrix related to the kinetic energy.On the right-hand side, q ∼ ext is the external acceleration flux imposed at the unit cell boundary, i.e. the acoustic load.The system (3.4) is symmetric.When a time-harmonic response is considered, it has real positive eigenvalues since no dissipative phenomena are included in the micro-scale constitutive models (3.2).The finite element nodes may be subdivided into the nodes at the boundary, where the essential boundary conditions are prescribed, denoted by 'ʙ', and the remaining 'independent' nodes that are internal to the unit cell, denoted by 'ɪ', resulting in the partitioned form of equation (3.4): Next, the goal is to express the internal degrees of freedom p ∼ɪ in terms of a reduced space.
For the considered wave problem, two distinct and approximately superimposable solutions can be discerned for constructing a representative reduced basis.The first solution assumes a steady-state response subjected to arbitrary Dirichlet boundary conditions on the ʙ boundary nodes: The solution to this problem is used to construct the so-called long wavelength basis, i.e. the unit cell response at the limit of a long wavelength propagation, enabling the elimination of the internal degrees of freedom p ∼ɪ .As a result, the solution of the transient unit cell problem (3.5) projected onto the long wavelength basis implies an instantaneous response of the internal pressure nodes p ∼ɪ according to the steady-state condition computed by a given t ∈ t 1 , t 2 .
The second solution is the eigenvalue problem where ω is the eigenfrequency in (rad s −1 ) of the unit cell subjected to zero Dirichlet boundary conditions.It gives the internal dynamics when the unit cell's outer boundary is close to a steady state relative to the internal domain, such as the case considered in this article.This condition was introduced in the context of computational homogenization of solid metamaterials [36,39,[45][46][47][48][49][50] and is called the relaxed separation of scales.The solution of the eigenvalue problem (3.7) provides the spectral basis, which using a modal expansion allows reducing the internal degrees of freedom p ∼ɪ to the ɴ lowest relevant eigenmode amplitudes η ∼ɴ , with typically ɴ≪ɪ.Indeed, the aim here is to capture the local resonance behaviour using a small spectral basis generated under zero Dirichlet boundary conditions.Sridhar et al. [36] demonstrated for a solid metamaterial case that the zero Dirichlet boundary condition leads to a suitable basis.Here, it will be shown that it is equally appropriate as a boundary condition for a fluid metamaterial described in the pressure formulation.Based on the above considerations, the solution to the problem of a fluid unit cell that has localized dynamics (subwavelength resonances) can be represented as the superposition of the long wavelength and local resonance solutions: (3.8) p where , m ʙɴ (c) and m ɴʙ (c) are the finite element submatrices of (3.5) projected on the reduced basis.Since the involved matrices are symmetric, it holds that m ʙɴ (c) = m ɴʙ (c) T .Note that the coupling between the local resonance dynamics (3.15) and the long wavelength dynamics (3.16), given by the matrices m ɴʙ (c) and m ʙɴ (c) , is purely driven by the inertial forces as these matrices depend only on the matrix Q and the reduced basis.

(c) The scale coupling
First-order spatial computational homogenization [37] expands the microscopic field variables in terms of the macroscopic fields and their gradients.For the acoustics problem considered here, the micro-scale pressure field p x → , t is represented as The micro-fluctuation field w p x → , t is added to accommodate the difference between the macro-scale field and the local variations that may develop within the unit cell due to the presence of inhomogeneities; x R → is a reference position vector, here chosen to be placed at the centre of the unit cell.
To couple the scales, in computational homogenization methods, an equivalence between the generalized power at the two scales is established, i.e. the so-called Hill-Mandel condition.To this end, the variational formulation of the macroscopic problem is extracted from equation (2.1).The integral functional is written as Using the divergence theorem, (3.19) can be rewritten as where n M → is the outward normal to the external boundary of the macroscopic domain Ω M .
In the left-hand side of equation (3.20), the first term is related to the elastic power rate, and the second term gives the kinetic power rate.The right-hand side of equation (3.20) describes the external virtual power rate applied to the macroscopic surface ∂Ω M .A similar procedure is performed for the micro-scale equation (3.1) on the unit cell domain Ω f .Then, the integral function of the micro-scale problem takes the same form as equation (3.20) with the additional Γ fs -internal interface term, on the left-hand side, which vanishes due to the interface condition (3.3), leading to where n f → is the outward normal to the outer boundary of the fluid unit cell domain Ω f .
Given the variational forms at macro-(3.20)and micro-scales (3.21), the link between the scales is recovered by equating the local macro-scale virtual power to the unit cell volume averaged power.Hence, the extended Hill-Mandel principle becomes The volume of the fluid domain V f is related to the total unit cell volume V by the porosity factor ϕ = V f /V.Using this definition and substituting (3.21) for the right-hand side of ( Substituting the macro-to-micro relation (3.18) into (3.23)gives As required in any computational homogenization approach, the virtual power of the microfluctuations on the boundary is forced to vanish by using appropriate boundary conditions, i.e. (3.25) For instance, this is ensured by zero uniform boundary condition on the micro-fluctuation, referred to as 'null boundary condition' by Blanco et al. [51] (3.26) Since equation (3.24) must hold for all admissible macro-and micro-fields, the macro-constitutive relations are obtained from the respective boundary integral relations, which can also be expressed as volume integrals by using the divergence theorem and the balance equation (3.1) For solid locally resonant metamaterials, the incorporation of the micro-inertia contribution into the macro-scale flux (stress) was introduced by Pham et al. [39].Here, an analogous expression is found for the macro-scale flux (acceleration) U → ¨M which depends not only on the micro-scale field U → ¨ but also on the micro-inertia term ϵ¨x → − x R → , see the right-hand volume integral of (3.27).This resulting expression is the transient version of the time-harmonic homogenization developed by Gao et al. [52] for a deformable porous solid saturated by a fluid.Finally, it is instructive to establish the scale coupling relations between the macro-and micro-scale pressure and gradient of pressure fields.The averaging relation for the pressure gradient can be found by applying the micro-scale gradient on (3.18) and integrating over the unit cell volume Ω f , followed by using the divergence theorem on the term involving micro- fluctuations.With regard to ansatz (3.26), the latter term vanishes, and the macro-to-micro constraint for the pressure gradient is derived Next, the averaging relation for the pressure can be found by integrating equation (3.18) over the unit cell boundary and taking x R → in the centroid of the unit cell, with account for equation (3.26), the terms ∂Ω f x → − x R → dA and ∂Ω f w p dA vanish.Then, the constraint for the pressure is obtained p M = 1 where A f = ∂Ω f dA.This macro-micro constraint aligns with the classic asymptotic and computational homogenization of locally resonant materials where the macro-scale field represents the averaged field of the non-resonant domain, e.g.see Boutin et al. [40], Pham et al. [39] and Sridhar et al. [36].

(d) Upscaling the reduced micro-scale response
Considering the infinitesimal acceleration flux dq ext = U → ¨⋅ n f → dA, for the discretized unit cell problem, the scale coupling relations (3.27) can be written as where 1 ∼ʙ is a column with scalar unit components.The notation Δ x η j ẅhere the upper-case homogenized material coefficients C ξ , D ξ , G ξ and H ξ , with subscript ξ = ϵ, U , correspond to the long wavelength regime, and the lower-case ones c j and d j → are the homogenized material coefficients of mode j, corresponding to the local resonance regime.Furthermore, combining equation (3.15) with (3.29) yields the evolution equation for the internal dynamic variables η N (3.31) where γ j and δ j → are the homogenized material coefficients of mode j, corresponding to the local resonance regime, where (3.32) γ j = Vc j , and δ j → = V d j → , for eigenmodes j = 1, 2, . . ., ɴ .
The detailed derivations leading to constitutive relations (3.30) and the evolution equation (3.31), as well as the final expressions of the homogenized coefficients, are given in appendix A, where it is also shown that the coefficients G ϵ , H ϵ → and G U → are identically zero.The long wavelength limits material properties C ϵ and H U are the conventional coefficients in the constitutive models in homogeneous (non-porous) fluid acoustics, and they are typically non-negligible, as expected since the homogeneous fluid is a limit case for ϕ 1.The effective macro-scale bulk modulus and mass density tensor are, respectively, Coefficients D U and D ϵ → = C U → are similar to the so-called Willis coefficients [24] and naturally 10 royalsocietypublishing.org/journal/rsta Phil.Trans.R. Soc.A 382: 20230368 emerge from the spatial unit cell averaging.The latter may or may not vanish depending on the unit cell symmetries and the specific choice of x R → , even in the long wavelength limit, which agrees with the findings of Willis [53].The local resonance coefficients c j and d j → depend, in general, on the nature of the occurring subwavelength resonance.For instance, a monopolar resonance induces a negative bulk modulus, leading to a large corresponding coefficient c j .On the other hand, if a dipolar resonance is present, it induces a negative effective mass, resulting in a large coefficient d j → .To summarize, at the macro level, an enriched homogenized continuum has emerged.The internal degrees of freedom p ∼ɪ , governed by equation (3.5), were projected onto the reduced space identified in §3b and upscaled as a superposition of the long wavelength and local resonance responses in §3d.Note that this model reduction and projection need to be done only once for a given unit cell, thus justifying the computational costs involved.The resulting formulation emerging at the macro-scale can be denoted as p M , η ∼ɴ , with p M and η ∼ɴ being the primary unknown fields.By virtue of the closed-form constitutive relations (3.30), the unit cell no longer has to be consulted, revealing a fundamental difference with respect to the FE 2 computational homogenization approach.To conclude, the enriched homogenized continuum, depicted in figure 1b, is governed by the conservation equation (2.1), the constitutive model (3.30) and the ɴ evolution equation (3.31) for η ∼ɴ enrichment variables.The latter can be interpreted to provide the enrichment of the homogenized continuum in a micromorphic sense.Figure 2 summarizes the multiscale procedure to construct the effective macro-scale model.

Numerical example
To illustrate the capabilities of the presented framework, an acoustic locally resonant metamaterial is considered, based on the working principles presented by Cheng et al. [8].The unit cell of size a contains a rigid solid labyrinthine-like structure, with size b, immersed in a fluid modelled in two dimensions, as shown in figure 3. The coiled space region triggers the fluid resonance inside the rigid structure in a subwavelength regime relative to the non-coiled region.For the unit cell in figure 3, when an acoustic wave in the non-coiled region travels a distance b/2, the wave inside the labyrinth travels a distance of approx.7b/2 through the labyrinth.Effectively, the wave speed is therefore reduced by a factor of 7. As discussed in §3a, given the dynamic viscosity of air ν f = 1.802 × 10 −5 m 2 s −1 and the smallest constriction d, the minimum frequency f min , for which the governing equation holds, is estimated to be 2 Hz. Figure 3b shows the two-dimen- sional finite element mesh used for the unit cell analysis.Quadratic serendipity elements are used with eight pressure nodes per two-dimensional element.

Air is considered as
Figure 4 shows the five lowest (local resonance) modes resulting from the eigenvalue problem (equation (3.7)) with zero Dirichlet boundary conditions.These modes constitute the spectral basis of size ɴ = 5.
Table 1 summarizes the constitutive parameters from (3.30) computed for the considered unit cell.Due to the symmetries of the unit cell, the homogenized fluid density ρ M and non- classic coefficient D U → are isotropic, while the Willis-like coefficients D ϵ → and C U → vanish.The presence of the rigid channels forces the fluid to move through a tortuous path, yielding an apparent increase of density, i.e. the homogenized fluid density ρ M is greater than the back- ground fluid density ρ f .On the other hand, for the homogenized bulk modulus K M , originating from the long wavelength limit basis, the equality K M = K f /ϕ holds, underlining the consistency of the framework with the analytically derived effective bulk modulus of a fluid with immersed rigid rods in the longwave limit as discussed by Hu et al. [30] and Li & Chan [3].
The dispersion curve of the enriched homogenized continuum can be calculated from the governing equations ( 2 p ^M The above system is Hermitian as follows from the relation between the off-diagonal terms given in (3.32), which is noticeable once the first line is pre-multiplied by the scalar constant −V (the unit cell volume).For a given real-valued wavevector k → spanning the contour of the irreducible Brillouin zone (M − Γ − X − M), the corresponding (real) eigenvalue ω 2 can be obtained by solving (4.1), which is plotted as a solid red line (-) in figure 5a.To validate this result, a direct numerical simulation of the fully resolved unit cell is conducted by applying Bloch-Floquet periodic boundary conditions, providing the reference solution for the dispersion spectrum, shown as blue crosses (+) in figure 5a.The Bloch analysis in this direct numerical simulation involves approx.ʙ+ɪ ≈ 6000 degrees of freedom.In contrast, the online stage of the homogenization approach contains only 6 degrees of freedom, corresponding to the size of the matrices in (4.1), which is 1 + ɴ = 6.It can be observed that the dispersion branch is very well predicted by the homogenized model.A discrepancy is observed in the second branch for frequencies above approx.550 Hz when the separation of scales is gradually violated.Frequency-dependent effective material properties can be obtained by assuming time-harmonicity and eliminating the enrichment variables η j in the constitutive model by inserting equation (3.31) into (3.30).Thus, the frequency domain effective homogenized constitutive relations read as Table 1.Homogenized material properties as defined in equation (3.30) computed for the unit cell of figure 3 (see also appendix A). local resonance local resonance where K M (eff) ω and ρ M (eff) ω are the homogenized effective bulk modulus and effective density, respectively.These are shown in figure 5b,c as a solid red line (-), whereas the long wavelength limit values, i.e. for ω 0, are shown as dashed black lines (−−).The first bandgap (red-shaded range around 400 Hz) is related to the effective bulk modulus, while the second bandgap exhibited by the homogenized model (red-shaded range around 700 Hz)   relates to a negative effective density.The effective Willis-like coupling coefficient S M → (eff) ω is approximately zero even in the vicinity of the resonances for the considered unit cell and is therefore not shown here.Interestingly, it is observed that the Willis-like coefficient has interaction term between monopolar and dipolar material properties.This is consistent with the findings of Sieck et al. [29] showing that the Willis coupling originates from the interaction between monopole and dipole motion from asymmetries at the micro-scale.The collection of effective material properties K M (eff) ω , ρ M (eff) ω and S M → (eff)  ω are analogous to the polarisability tensor of an individual scatterer described in [54][55][56].Appendix B shows the Willis coupling coefficient of an asymmetric unit cell to illustrate the framework's capabilities.
Next, the homogenized continuum's response in a finite-sized domain is evaluated on the configuration depicted in figure 6.A plane wave is excited at the left of the air acoustic domain.The wave enters the metamaterial layer composed of 5 unit cells sequentially arranged along the horizontal direction.The infinite character in the vertical direction is achieved by defining periodic boundary conditions on the top and bottom surfaces.Perfectly matched layer (PML) domains are employed to simulate infinitely long incoming and outgoing acoustic domains in the horizontal direction.The reference direct numerical simulation solves the full-field acoustic problem, as shown in figure 6a.
The homogenized response is computed using the model illustrated in figure 6b, combining a fluid (acoustic) domain in a p a -formulation in white colour coupled to the homogenized enriched continuum domain in a p M , η ∼ɴ -formulation in grey colour, implemented in COMSOL ), obtained via the computational homogenization approach given by the expressions defined in equation (4.2).
where n → represents the normal to the interface between the enriched continuum and the plain air.The mesh of the enriched continuum, depicted in figure 6b, consists of quadratic serendipity elements with a size equal to the unit cell size, i.e. 5 elements.The number of degrees of freedom for the full-field model depicted in figure 6a was approx.40 000, while the model including the homogenized continuum of figure 6b had around 350 degrees of freedom, representing a reduction of approx.99.1%.To reduce computational costs, time-harmonicity of the primary fields is assumed, i.e. the simulations are performed in the frequency domain (steady-state).However, it should be stressed that the proposed framework-the initial boundary value problem (2.1), (3.30) and (3.31)-can equally be used to solve transient problems on arbitrary geometry in the time domain, i.e. problems that are not necessarily a three-medium transmission as considered here for illustration purposes.Following the approach by Lewińska et al. [50,58], a three-microphone method is employed to compute the transmission loss through the metamaterial layer.For a given excitation frequency ω, the complex pressure amplitudes p 1 , p 2 and p 3 are obtained as the pressure p a ( x → , ω) averaged along the vertical lines at positions x 1 , x 2 and x 3 , as indicated in figure 6a,b .The reflection and transmission coefficients are calculated, and the transmission loss is determined using the expressions: Here, L = 5a is the total layer thickness, and k a = ω/ K a /ρ a is the wavenumber in the horizontal direction in the plain air acoustic domain.The transmission losses, computed with the direct numerical simulation and the homogenized enriched continuum, are compared in figure 6c.An excellent agreement is observed up to approx.550 Hz.However, a discrepancy between the homogenized continuum and the reference solution is evident for higher frequencies from approx.350 Hz onwards, as also observed in the dispersion analysis (figure 5).The underestimation of the sound speed shown in figure 5 correlates with the underestimation of the transmission loss peak due to resonance in figure 6.This can be attributed to the decomposition (3.8), which is no longer valid for a short wavelength with respect to the unit cell size, including the largest peak due to the negative bulk modulus effect.A study on the influence of the macro-scale mesh has been performed (not shown here).It has been found that for a metamaterial with a macro-scale mesh with more than five elements, the transmission curve remains practically unchanged.On the other hand, a significant error is observed for fewer than five elements, which is to be expected since, in this case, the effective wave is not accurately discretized.
As evident from the above results, the framework is valid for (low) frequencies-up to approx.550 Hz in the considered case-that satisfy the scale separation at the unit cell boundary, i.e. the non-coiled region, which can be referred to as the host domain.Therefore, the effectiveness of the framework for accurately upscaling the micro-scale dynamics relies on the subwavelength character of the resonance, which is given by the contrast between the host domain (non-coiled) and the resonant domain (the labyrinth).For the considered unit cell, there is a geometric contrast of approx.7 times, whereas, for instance, in the unit cell of Liu et al. [1], there is a material contrast of 1000 in stiffness between the host matrix and the resonant rubber-coated inclusion, ensuring a highly localized (subwavelength) resonance.To conclude, high contrast between host and resonant domains leads to a unit cell with more locally resonant modes that can be upscaled with the present framework.Here, a relatively low contrast is analysed, and it illustrates well the limit of the hypothesis on the validity of the equation (3.8), which is satisfied only by the first mode in this case.

Conclusion
The high computational cost for finite-size metamaterial simulations arises from its fine geometric features.Reduced-order computational homogenization is a solution to address this issue.This article proposes an efficient multiscale model for analysing the transient response of acoustic LRAMs based on a pressure formulation and it is applicable to both frequencyand time-domain analysis.The method couples macro-and micro-scales using an extended version of the Hill-Mandel principle.The micro-scale dynamics are represented in a reduced space as the superposition of long wavelength and local resonance responses.The outcome is a macroscopic homogenized enriched continuum with additional enrichment variables describing the micro-scale dynamics in a subwavelength regime.The enrichment variables may be eliminated if a frequency domain analysis is performed.The upscaled effective constitutive model is versatile, modelling exotic macroscopic responses typically found in acoustic metamaterials.The dispersion relation and the considered finite-size problem indicate that the homogenized model captures the behaviour of the direct numerical simulations very well within the limits of its applicability, demonstrating the framework's efficacy and accuracy.
Other versatile results can be obtained with the presented framework.For example, a thorough investigation of the unit cell symmetries and their consequences on the homogenized coefficients for an LRAM unit cell with non-negligible Willis coupling.Another result is a straightforward extension of the current framework considering a deformable solid internal to the unit cell, i.e. disconnected solid structures such as the unit cell of this article.A solid and fluid resonance can be designed to take place in a close frequency range, yielding interesting metamaterial behaviour originating from co-existing subwavelength resonances.Notably, in 16 royalsocietypublishing.org/journal/rstaPhil.Trans.R. Soc.A 382: 20230368 the latter case, the macro-scale model remains the same as presented here: an equivalent fluid supporting only compression waves since the solid phase is not connected.
The approach developed in this article will facilitate the design of metamaterials in realistic industrial scenarios where fast simulations are required for optimization.Therefore, the proposed homogenization framework is an important step to accelerate progress towards achieving a higher Technology Readiness Level for the emerging class of metamaterials, where numerous applications await exploration of their full design space, as pointed out by Kochmann et al. [59].Moreover, this framework paves the way for developing multiscale methods for high-fidelity acoustic problems, where more elaborate material models are required, including dissipation and flow-induced nonlinearities; for instance, acoustic liner applications, as discussed in Jones et al. [60] and Rego et al. [61].
The coefficients G ϵ , H ϵ → and G U → depend on the product H ʙʙ 1 ʙ , vanishing identically, since 1 ʙ belongs to the null space of H ʙʙ .Mathematically analogous to the solid rigid body motion, this corresponds to the case when the fluid undergoes a homogeneous change in pressure 1 ʙ on the prescribed nodes ʙ, in the absence of flow, in other words, a homogeneous volume expansion/contraction takes place.Neglecting the identically zero coefficients, the eigenvalue problem reads ).It, therefore, preserves the self-adjointness of the operator on the state vector ^M η ∼ɴ T that leads to an eigenproblem with real eigenvalues, which is expected for a conservative system.Logically, the eigenvalue problem of the homogenized continuum has real eigenfrequencies since the full-field response does not decay with time, which aligns with the results of Willis [63] showing that self-adjointness at the micro-scale implies self-adjointness at the macro-scale.
(a) Governing equations Consider an acoustic metamaterial unit cell where a fluid is constrained to move through a rigid solid structure, as depicted in figure 1c.The acoustic wave equation governs the fluid dynamics 4 royalsocietypublishing.org/journal/rsta Phil.Trans.R. Soc.A 382: 20230368 linearized around a thermodynamic equilibrium state, assuming the fluid is an ideal gas with a quiescent background, (3.1) −ϵ¨x

Figure 1 .
Figure 1.(a) Full-field problem.A rigid solid domain (black) and a connected compressible ideal fluid domain (white) define an acoustic metamaterial domain subjected to boundary constraints and time-dependent acoustic loading.(b) Equivalent macroscopic domain Ω M subjected to boundary constraints and Dirichlet time-dependent acoustic loading p D t with a representative macroscopic wavelength λ M .(c) A unit cell with the definition of the micro-scale rigid solid (black) domain Ω s enclosed by fluid-solid interface Γ fs .The fluid (white) domain Ω f is enclosed by the surface ∂Ω f and fluid-solid interface Γ fs .The metamaterial unit cell may have resonant cavities such as the depicted Helmholtz resonator.

j 2 ,
/journal/rsta Phil.Trans.R. Soc.A 382: 20230368The discrete form of the long wavelength limit equation(3.6) is obtained by disregarding the transient term in equation(3.5).From there, p ∼ɪ (lw) can be expressed in terms of the long wavelength basis Σ ɪʙ as (3.9) p ∼ɪ(lw) = Σ ɪʙ p ∼ʙ with Σ ɪʙ = − H ɪɪ −1 H ɪʙ .The discretized time-harmonic problem (3.7) is given by (3.10) −ω j 2 Q ɪɪ + H ɪɪ φ the jth eigenmode, where j = 1, …, ɪ.Next, a selection of ɴ lowest localized resonant modes will be made, with the truncation ɴ≪ɪ.The selected modes are collected in the matrix Φ ɪɴ referred to as the 'local resonance' basis.Thus, the internal nodal pressure variables in the local resonance solution are approximated using the truncated spectral basis (3.11) p ∼ɪ (lr) ≈ Φ ɪɴ η ∼ɴ with η ∼ɴ the column of modal amplitudes.The matrix Φ ɪɴ is normalized with respect to the volumetric strain matrix, i.e. satisfying (3.12) Φ ɪɴ T Q ɪɪ Φ ɪɴ = I ɴɴ (lr) and hence Λ ɴɴ (lr) = Φ ɪɴ T H ɪɪ Φ ɪɴ with I ɴɴ (lr) the identity matrix of the appropriate dimension.As a direct consequence, the second equation of (3.12) is the kinetic energy matrix expressed in the modal basis.The diagonal matrix contains the truncated eigenvalues ω j where j = 1, …, ɴ.Collecting (3.9) and (3.11) into (3.8), the transformation A resulting from the superposition of the solutions projected on the reduced bases I are the zero and identity matrices, respectively.Equation (3.5) can then be projected on the reduced basis(3.14) the background fluid with ρ f = 1.225 kg m − 3 and K f = 1.441 × 10 5 Pa at a temperature of 15 ∘ C. The values of the geometric features indicated in figure 3a are a = 105, b = 72.7,d = 4.05, ℎ = 31.3and w = 1 in mm.

∇Figure 2 .
Figure 2. The summary of the reduced-order computational homogenization steps for constructing an equivalent macro-scale model from the selection of a unit cell.

Figure 3 .
Figure 3. (a) Locally resonant labyrinthine metamaterial with a rigid solid structure that coils the wave path.(b) Finite element mesh for the fluid domain.

Figure 4 .
Figure 4.The lowest five eigenmodes constitute the truncated spectral basis with ɴ = 5.The coloured domain is the fluid phase, whereas the rigid solid domain is shown in black.

Figure 5 .
Figure 5. (a) Dispersion curves obtained through the computational homogenization approach plotted as a solid red line (-) versus the dispersion curves computed using the direct numerical simulation Bloch analysis (DNS) plotted as blue crosses (+).(b) Effective bulk modulus and (c) effective mass density, component e 1 → ⊗ e 1 → (equals to component e 2 → ⊗ e 2 →), obtained via the computational homogenization approach given by the expressions defined in equation (4.2).

14 royalsocietypublishing
.org/journal/rsta Phil.Trans.R. Soc.A 382: 20230368 using the weak form PDE interface in the Mathematics module.The subscript a denotes the plain air domain, which is not homogenized and has material properties ρ a = ρ f and K a = K f .The coupling between the enriched continuum and the plain air domain is established by imposing continuity of pressure and continuity of normal acceleration, as discussed in Allard & Atalla[57] for porous materials (4.3) p M = p a and ϕU

Figure 6 .
Figure 6.(a) Full-field model constituted of rigid solid labyrinthine structures (black) immersed in a plain air domain (white).(b) The equivalent macroscopic model where the metamaterial layer is considered as a homogenized continuum (grey) coupled to a plain air domain (white).The macro-scale mesh is shown only in the homogenized domain.(c) The transmission loss curve was computed via the three-microphone method for models (a) and (b).
the choice of x R → is made such that certain coefficients are significantly reduced, for example, D U and (D ϵ → , C U → ), in which the latter pair becomes zero in the case of a square symmetric unit cell.Note that even in the case (D ϵ → , C U → ) are nonzero, for an asymmetric unit cell, the sum −D ϵ → ⋅ k → + k → ⋅ C U → vanishes in the diagonal of matrix (A 5).This follows from the equality of the coefficients D ϵ → = C U → since Q ʙʙ is a symmetric matrix (see equations (A 3) and (A 4)

Figure 7 .
Figure 7.The absolute value of the effective Willis coupling coefficient obtained via the computational homogenization approach for a symmetric unit cell plotted as a solid red line and an asymmetric unit cell plotted as a solid blue line.The long wavelength limits are dashed black lines.Interpretation of curve's colors in the figure is not necessary since the corresponding unit cell is indicated with an arrow.