Computational homogenisation of acoustic metafoams

Acoustic metafoams are novel materials recently proposed for low frequency sound attenuation. The design of their microstructure is based on the combination of standard acoustic foams with locally resonant acoustic metamaterials. This results in improved sound attenuation properties due to the interaction between viscothermal dissipation effects and the local resonance effects at the pore level. In this paper, the non-standard behaviour of such a metafoam with a complex two-phase microstructure is analysed through a multiscale approach. The macroscopic problem is described by general balance equations and at the microscopic scale a detailed representation of the microstructure is considered. The frequency dependent effective properties are used to explain the extraordinary acoustic performance. The homogenisation approach is also validated using direct numerical simulations, showing that the homogenisation technique is adequate in modelling both viscothermal dissipation and the local resonance effect within the metafoam microstructure.


Introduction
Low frequency sound attenuation in a subwavelength regime is a challenging task. On one hand, the mass law dictates a heavy weight for a sound isolation barrier (Ballou, 2015), whereas on the other hand, in order to perturb or absorb the waves of long wavelengths, a large thickness of such a barrier is required (Long, 2005). Recently, in the field of acoustic metamaterials, materials and structures are designed with exotic properties, uncommon or non-existent in nature, overcoming constraints for sound attenuating materials. Among potential solutions for low-frequency shielding/attenuating purposes, several approaches have been proposed: fractal acoustical metamaterials (Song et al., 2016) with improved absorptive and reflective properties due to the presence of multiple resonances, dark acoustic metamaterials (Mei et al., 2012) characterised by broad-band total absorption obtained with an elastic membrane decorated with asymmetric rigid platelets or ultra-sparse metasurfaces (Cheng et al., 2015) providing a high reflectance due to either a negative bulk modulus or a negative mass density, to name a few.
A different approach towards the low frequency noise challenge has recently been presented in Lewińska et al. (2019), where instead of designing a metastructure with an extraordinary behaviour, a new class of materials, acoustic metafoams, has been proposed. Acoustic metafoams rely on the combination and interaction of two interplaying attenuation mechanisms, namely, the viscothermal dissipation occurring in standard acoustic foams and local resonance effects. This combination results in enhanced wave absorption and reflection. The demonstration of this concept was based on a poro-elastic unit cell with an embedded resonating mass, where the numerical studies involved band structure analyses and transmission spectrum calculations. However, due to the geometrical complexity of the microstructure and its small size, the direct numerical simulations of transmission were restricted to a thin material layer (up to 15 unit cells). Therefore, there is a strong need for an efficient computational approach to assess the performance of such materials in realistic finite size applications.
Homogenisation techniques are natural candidates for improving the computational efficiency, allowing at the same time to account for the microscale morphology as well as the complex physical phenomena occurring at the microstrucural level. Typically, the relationship between the microstructure of a composite and the macroscopic performance of the material is established through frequency dependent effective parameters. Within a dynamic homogenisation framework, the exotic behaviour of acoustic/elastic metamaterials manifests itself in a non-standard frequency dependence of the effective mass density and elastic moduli, including the variation of these parameters towards negative values (Brunet et al., 2013). Starting from a discrete description of locally resonant acoustic metamaterials, an analytical model can be used to obtain the frequency dependent mass density in 1D (Liu et al., 2005) and in 2D (Huang et al., 2009), and the effective medium theory can be used to predict double negative behaviour (Wu et al., 2007). In the literature, some useful developments for the modelling of wave dispersion and dissipation in heterogeneous periodic materials https://doi.org/10.1016/j.euromechsol.2019.103805 Received 2 May 2019; Received in revised form 11 June 2019; Accepted 13 June 2019 have been proposed. Just to mention a few, Bloch based methods proposed by Willis (2012) have been used to solve inhomogeneous elastodynamics problems in steady state (Nassar et al., 2015). Among the recent contributions within enriched homogenisation schemes, the computational homogenisation framework proposed by Sridhar et al. (2016Sridhar et al. ( , 2018 and the variational homogenisation framework developed by Roca et al. (2018) are to be mentioned, since both can be adopted for transient modelling of locally resonant acoustic metamaterials (both including micro-inertia effects).
The behaviour of foam-like materials is commonly described using Biot's theory (Biot, 1956). However, for poro-elastic materials, for which the skeleton motion affects the fluid flow within the microstructure or the microscopic solid deformation is significant, Biot's approach is no longer sufficient. In such cases, asymptotic homogenisation frameworks have been adopted for instance by Yamamoto et al. (2011) for regular poro-elastic microstructures or by Venegas and Boutin (2017) for permeo-elastic materials.
In turn, Gao et al. (2016) proposed a computational homogenisation method based on general balance laws which constitute an extension of Biot's theory. The method of Gao et al. is a multi-scale approach, for which the macroscale problem is described by the classical equations of momentum and mass balance, whereas at the microscale, instead of using the rule of mixtures (Schraad and Harlow, 2006;Coussy et al., 1998) or ensemble averaging (Lafarge and Nemati, 2013), a detailed description of the unit cell with the complete fluid and solid domains is provided through which the visco-thermal effects are accounted for. Recently, this homogenisation framework has been successfully used to investigate the influence of the microstructure on the performance of acoustic foams (Gao et al., 2017).
The aim of this work is to investigate the wave attenuation behaviour of acoustic metafoams through a computational homogenisation approach. To this end, the framework proposed by Gao et al. (2016) is used as a point of departure. The main novel aspects of this work are: • identification of the frequency dependent effective material parameters of acoustic metafoams and revealing the relation between microstructural properties and attenuation performance of the material; • analysis of emerging effects due to the interaction between the local resonance and the fluid viscosity with respect to the effective parameters and acoustic behaviour of metafoams; • demonstrating the applicability of the homogenisation framework for modelling acoustic metafoams and predicting sound attenuation effects for engineering applications.
The paper is organised as follows. In Section 2, the computational homogenisation scheme is presented followed by its particularization for to the acoustic metafoams in a transmission set-up. In Section 3, results are shown with focus on a comparison between direct numerical simulations and homogenised models for several unit cell configurations including or excluding local resonance effects and with different fluid viscosities. The paper ends with concluding remarks.

Computational homogenisation for poro-elastic media
In this section, the main points of the multiscale computational homogenisation framework proposed in Gao et al. (2015Gao et al. ( , 2016 are presented. Two coupled problems are considered: the macroscopic problem where the poro-elastic material is replaced by a homogeneous material and the microscopic one with a fully heterogeneous representative volume element (RVE) consisting of fluid (air) and solid (polymer, metal) domains. It is assumed that scale separation holds, meaning that the characteristic length of the external excitation is much larger than the microscopic characteristic length (i.e. the wavelength is much larger than the cell size). The studies at both levels are conducted in the frequency domain adopting the + i convention with i denoting the imaginary unit and ω being the angular frequency.

Microscopic problem
The microscopic poro-elastic RVE consists of a solid and a fluid domain. For the sake of simplicity, no specific subscripts are used to indicate the variables associated with this scale. The solid constituent is modelled as a linear elastic material with the momentum balance equation given by: where u s is the solid displacement vector, s 0 the solid density and s is the stress tensor given by = C u : s s 4 , with C 4 being the fourth-order elasticity tensor, which is assumed isotropic. Symbol denotes the gradient operator at the microscopic scale. Note, that unlike in Gao et al. (2016), here, the thermal effects are not included in the description of the solid domain which is justified by the high contrast in thermal conductivity between solid and fluid. The solid is considered under isothermal conditions at the ambient temperature.
The fluid domain is described using the linearised Navier-Stokes-Fourier equations: which represent the balance of momentum, the energy balance and mass conservation (considering the ideal gas law), respectively. In these equations, v f denotes the velocity vector which can be expressed in terms of the displacement as: 0 the equilibrium density, c p f the heat capacity at constant pressure and P 0 and 0 the ambient pressure and temperature, respectively.
Constitutive equations which complete the description of the fluid domain are: where µ f is the viscosity, f the thermal conductivity and I the second order identity tensor. The coupling between fluid and solid domains at the interface S i is prescribed by continuity of velocities and tractions (with n i being a normal vector pointing from the fluid to the solid domain): along with an isothermal boundary condition: At the outer boundaries of the RVE, a periodic boundary condition for the solid displacement is applied on the solid surface S e s : and a prescribed traction is defined on the fluid surface S e f :

Macroscopic problem
The macroscopic problem is described in terms of the solid displacement u M s and the fluid pressure p M f . The macroscopic fluid viscous stress is assumed to be negligible compared to the macroscopic pressure and sound propagation is considered as an adiabatic process. The porosity ϕ is defined as the volume fraction of the fluid and remains constant. Therefore the governing equations, momentum conservation for the solid and mass conservation for the fluid, are given by: where M s is the macroscopic stress tensor in the solid, f M s is the inertial force exerted on the solid,  (12) and (13) are solved together with macroscopic constitutive equations and appropriate boundary conditions.
Following Gao et al. (2015Gao et al. ( , 2016, a set of constitutive equations derived from the solution of the microscopic problem through homogenisation is given by: where the fourth order tensor C 4 is the stiffness tensor dominated by the stiffness of the solid skeleton and the scalar S f is dominated by the compliance of the air. The symmetric second-order tensor sf is the coupling parameter and M K K , , s f sf are also second-order tensors. The effective properties in equations (14)-(17) are required to solve the macroscopic problem and are obtained through a number of microscopic simulations.
At the macroscopic level a relationship between the homogenisation framework proposed by Gao et al. and Biot's theory can be established. For this purpose, Eqs. (16) and (17) can be rewritten as: (20) where av s denotes the averaged solid density (for a possibly inhomogeneous solid phase). Assuming that the microscopic solid displacement fluctuations can be ignored and the volume-average fluid deformation is small enough, the following approximations hold: Under these assumptions the homogenised equations reduce to Biot's theory.

Micro-to-macro relations
The relationship between micro and macroscale is established through the Hill-Mandel condition, i.e. the energy variation per unit volume in a macroscopic point equals the volume average of the total microscopic energy variation per unit volume of the associated unit cell: The macroscopic virtual energy rate can be written as: and the volume average of the microscopic virtual energy rate is: where n e s and n e f are the outward normal vectors on the solid and fluid surfaces, respectively. Recalling the energy consistency principle (25), the following micro-to-macro relations result: = By applying Gauss's theorem and using the balance equations of the microscopic problem, the boundary integral form of equations (28)-(31) can be rewritten by means of volume integrals: Note that equation (32) shows that both fluid and solid microscopic inertial terms participate in the macroscopic solid inertia, whereas equation (33) demonstrates that the macroscopic solid stress depends on the microscopic stresses in both domains and on both microscopic inertial contributions, which also introduces an intrinsic length scale in the macroscopic problem, related to the unit cell size, as can be observed from the factor x x ( ) R . This aspect has been pointed out previously, where it was discussed in the context of the transient computational homogenisation scheme by Pham et al. (2013). Equations (34) and (35) show that microscopic volumetric changes of both the solid and fluid phase influence two macroscopic quantities, i.e. the volumetric fluid change and fluid displacement.

Unit cell description
In order to assess the applicability of the homogenisation framework to an acoustic metafoam, a cubic unit cell is considered, following Lewińska et al. (2019). This unit cell is a simplified representation of a single foam pore along with an embedded resonator, constituting the basic unit cell of the acoustic metafoam.
The unit cell consists of a solid polyurethane (PU) skeleton made out of thick struts in a cubic arrangement. The struts are connected by thin membranes among which four out of six are partially open. The rest of the unit cell is filled with fluid (air). A resonator is additionally attached to one of the membranes in the form of a cantilever beam. Based on Lewińska et al. (2019), two unit cell configurations are considered. The first one, denoted as "concentrated mass", is characterised by a heavy mass at the tip of the resonator. The second configuration, denoted as "without mass", is characterised by the presence of the resonator without any added mass. These two configurations are shown in Fig. 1a and b. In addition to the configurations analysed in Lewińska et al. (2019), a third configuration (denoted as "distributed mass") it is considered here, which is a unit cell characterised by a homogeneously distributed mass over the solid domain equivalent to that of the unit cell with extra mass added. This last configuration is used to assess the influence of mass distribution on the performance of the unit cell.
The following dimensions for the pore are assumed. The size of the unit cell is = a 100 μm, with a strut thickness of 25 μm and a membrane thickness of 1 μm. The membrane opening is assumed to be 1% of the membrane surface. In the case of the unit cell with the added concentrated mass, a mass of 10 10 kg and dimensions 10 μm×10 μm×10 μm is attached to the tip of the cantilever resonator. A cantilever mass is used in order not to enhance the dissipation by an increase of the solid fluid interface. The volume of a high density metal material corresponding to this mass would not exceed the volume of the cell. Material properties for both the solid and fluid domains are listed in Tables 1 and 2. In the case of the distributed mass, the solid density is = av 1558kg/m 3 with no changes in its elastic constants (Table 2).
In the unit cell simulations, a periodic geometry is adopted. The boundaries of the domain used in the simulations are taken at a distance of a 0.2 relative to the position of the membranes, as illustrated in Fig. 1b and c. With this choice, the regions with high gradients of pressure and velocity (the membrane opening and the resonator itself) are located inside the unit cell, distant from the boundaries. As a consequence, for such a unit cell, the viscous stress contribution at the fluid outer boundary can be neglected as done by adopting Eq. (11).

Calculation of effective properties
In order to obtain the effective properties for the constitutive equations (14)-(17) for each of the considered unit cells, nine simulations with different loading conditions are conducted, owing to the symmetries of the unit cells (symmetry in the xy plane, making the x and y directions identical, but z different). The specific loading conditions for each simulation set are listed in Tables 3 and 4. The exact arrangement of the effective matrices is given in Appendix.
In order to assess the influence of the resonating mass on the behaviour of the material, the effective properties obtained for the unit cell with the concentrated mass are compared to those obtained for the cases without extra mass and with the distributed mass. The components of the tensors C 4 , sf and the scalar S f for these three microstructures have small imaginary parts, which are equal and almost constant in the frequency range considered. This means that the mass distribution within the unit cell does not influence the effective stiffness of the unit cell and the dissipative effects captured by these set of parameters are low. For instance, in Fig. 2, the effective dynamic fluid bulk modulus = B S / f f normalised by the ambient pressure P 0 is presented. It can be observed that the fluid compliance is constant in frequency (the real part) and thermal dissipation effects are low (the imaginary part). The latter remark is coherent with the findings in Lewińska et al. (2019), where the dissipative mechanism for this type of microstructure was shown to be dominated by the fluid viscosity.
The tensors M s and K sf (where (M s 1 , M s 2 , M s 3 ) and (K sf 1 , K sf 2 , K sf 3 )  Table 1 Material properties of air ( f 0 denotes fluid equilibrium density, f thermal conductivity, c p f heat capacity at constant pressure, µ f dynamic viscosity, 0 and P 0 equilibrium temperature and pressure, respectively) (Deckers et al., 2015).
1.2 0.0257 1005 1.84 10 5 293 1.01 10 5 Table 2 Material properties of the polyurethane (PU) solid phase ( s denotes solid density, E Young's modulus, ν Poisson's ratio) (Gao et al., 2017). Lewińska, et al. European Journal of Mechanics / A Solids 77 (2019) 103805 denote the diagonal components) are the effective properties that are influenced significantly by the presence of the resonator with the concentrated mass. In Fig. 3a, the first component M s 1 is shown, normalised with respect to the static densities of = av 277 kg/m 3 (volume-averaged over the entire unit cell) for concentrated and distributed cases and = av 177 kg/m 3 for the case without extra mass. The parameter M s 1 can be interpreted as an effective dynamic solid density, which exhibits a large positive maximum and large negative minimum around the resonance frequency.
Characteristic for locally resonant acoustic metamaterials, within the (shaded) band gap regime estimated in the range 445-560 Hz by Lewińska et al. (2019), the effective mass density becomes negative. At the resonance frequency of 445 Hz, also a negative dip in the imaginary part of M s 1 can be observed. Its sharp character is the consequence of a fully elastic response of the solid material (Park et al., 2012) as well as the moderate value of the fluid (air) viscosity. Note, that in the long wavelength limit, the values of the effective mass density converge to its static (volume-averaged) quantities. For the two other cases (without the extra mass and with the distributed mass), the effective densities M s 1 do not depend on frequency and remain equal to the volume-averaged values.
The behaviour of the parameter K sf 1 depicted in Fig. 3b for the three cases follows the same trend, i.e. it peaks near the resonance frequency in its real part and reveals a negative dip for its imaginary part for the unit cell with the concentrated mass. The parameter K sf 1 stays purely real and shows only a weak frequency dependence for the other two cases (a slight growth with frequency can be noticed).
In Fig. 4, the frequency variation of the parameter K f 1 , which is related to the fluid permeability, is depicted. Both real and imaginary values of this parameter are almost identical for all three cases. For the cell with a concentrated mass, around 445 Hz, the frequency dependence of Im K ( ) f 1 slightly deviates from being linear; however, this deviation does not exceed 0.2 10 5 .
The impact of the fluid viscosity on the effective properties M s 1 , K sf 1 and K f 1 obtained for the unit cell with concentrated mass is depicted in Figs. 5 and 6. With decreasing viscosity, the peaks and dips of the real and imaginary spectra of both M s 1 and K sf 1 become sharper. In turn, the higher the viscosity, the less pronounced is the parameter variation within the band gap region. The viscosity also affects the spectrum of parameter K f 1 . In particular, for the lowest considered viscosity = µ µ 0.1 air , the real part of K f 1 becomes closer to 0 and a significant increase of the imaginary part can be noticed (Fig. 6).
In order to verify if Biot's theory can be applied to describe the behaviour of an acoustic metafoam, the parameters s and f are calculated using Eqs. (20) and (21) and compared with the approximation given by Eqs. (23) and (24). In Fig. 7, the real and imaginary parts of parameters /(1 ) s 1 and / f 1 (which are the first diagonal components) are shown for the cases with the distributed and concentrated mass. Clearly, the approximation underlying Biot's theory does not apply in both cases (within the entire frequency range considered, the real values are not equal to 1). Note that for the concentrated mass case, the peak and dip around 445 Hz are still present, which is a consequence of the localised resonance, as discussed before.

Macroscopic simulations
The homogenisation model is evaluated next through the comparison with direct numerical simulation (DNS) results of a transmission test (Lewińska et al., 2019). The transmission set-ups shown in Fig. 8 mimic the experimental transmission impedance tube test used for assessing the acoustic performance.
In both cases, the sample is embedded in a homogeneous acoustic domain, described by the Helmholtz equation with properties of air: wave speed = c 343 m/s, density = 1.2 kg/m 3 . Plane wave excitation with unit amplitude is applied on the left boundary, whereas on the right, a perfectly matched layer is used in order to reduce spurious reflections. On the lateral planes, periodic boundary conditions are assumed and interface conditions between the acoustic domain and the homogenised metafoam layer are applied following Debergue et al. (1999). A three point method (Ho et al., 2005;Henríquez et al., 2017) is adopted to obtain transmission, reflection and absorption spectra based on the averaged pressure values evaluated at the three positions indicated by the red planes in Fig. 8.
In Fig. 9, transmission, reflection and absorption plots are shown for both the homogenised and DNS models for the three unit cell cases considered. First of all, it can be observed that the attenuation in the distributed mass case is in general stronger than in the case without an extra mass due to the mass law and it occurs mostly through reflection. Secondly, comparing the cases with the same total mass, it can be seen that the concentrated mass case exhibits a higher attenuation at frequencies around 455 Hz (a transmission dip occurs) than its distributed counterpart. The reflection peak in the concentrated mass case is however followed by a reflection dip, which leads to mitigation of the Table 3 Loading conditions used to determine the tensors M K K , , s f sf , with vectors e x , e z denoting unit vectors in x and z directions, respectively (see Fig. 1 Table 4 Loading conditions used to determine parameters S C, ,     attenuation performance at higher frequencies (around 600 Hz), which is typical for resonance based materials. The absorption obtained for the unit cells with a concentrated mass is locally enhanced in comparison with the other two cases.
In general, an adequate correspondence between the homogenisation and DNS results can be observed for all three cases. Considering the homogenised model and DNS results obtained for the unit cell with a concentrated mass, an almost perfect match can be observed, in particular, at frequencies in the range of 400-600 Hz. The homogenisation framework is therefore able to capture the local resonance effect and its influence on the acoustic performance of the material since both reflection and absorption spectra are in line with the DNS predictions. The ability of the homogenisation framework to reproduce this behaviour is due to the frequency dependent effective properties encoded in the constitutive equations. The presence of the reflection peak is a consequence of the peak and dip in the spectrum of the real part of the solid effective density M s 1 , along with very low values of the parameter Re K ( ) f 1 . On the other hand, the absorption level is governed by the imaginary parts of the parameters M K K , , s sf f 1 1 1 . A minor discrepancy between the homogenised and DNS models exists for all cases, which increases with frequency. This is related to the gradual deterioration of the assumed scale separation assumption; in other words, the higher the frequency, the more important secondorder effects will become. Moreover, for higher frequencies, the viscous stress contribution becomes more significant, whereas the boundary conditions adopted for the RVE as given by Eq. (11) neglect the viscous stresses. In Fig. 10, the ratios between the effective viscous stress and  the local pressure averaged on the cross-sections in the DNS for the case without extra mass are shown. Ten cross sections for each shift are considered, where the shift 0 denotes cross sections through the cell membrane and the shift a 0.5 is related with cross sections taken through the resonator. Based on Fig. 10, it can be stated that for the adopted choice of RVE periodicity (with the shift of a 0.2 ), the viscous stress to pressure ratio is relatively low, however, increasing with frequency. A similar conclusion holds for the other unit cell cases considered.
As mentioned in Lewińska et al. (2019), the wave attenuation in case of acoustic metafoams emerges from the presence of the local resonance together with the strong coupling between solid and fluid domains. A variation of the fluid viscosity (which determines the level of coupling) significantly influences the effective parameters within the resonance frequency range (for the concentrated mass unit cell), and as a consequence, the homogenised acoustic performance at the macroscale is affected as shown in Fig. 11. The highest attenuation at a frequency of 445 Hz occurs for the (default) air viscosity, however the broadest attenuation range is exhibited for a higher viscosity = µ µ 10 air . In fact, for the higher and lower viscosity cases, the attenuation mechanism changes relative to the air viscosity case; the reflection peak is reduced and higher levels of absorption emerge. The origin of this macroscopic behaviour is embedded in the effective properties. In the higher viscosity case, the lower reflection level is related with a significantly smoother spectrum of the solid effective density M s 1 . The absorption peak emerges as a result from the more persistent smooth imaginary part of parameter K sf . The weaker reflection for the lower viscosity case is caused by a higher value of Re K ( ) f 1 and the higher absorption is the effect of an increase of Im K ( ) f 1 . These findings can be directly translated to the unit cell design, where instead of modifying the viscosity of the fluid, analogous effects can be achieved by for instance changing the membrane opening ratios.

Conclusions
In this paper, for the first time a computational homogenisation approach has been applied to study the behaviour of acoustic metafoams. These materials are challenging from the homogenisation perspective, since not only the microstructure consists of fluid and solid parts with a strong coupling between them, but also a local resonance of the solid occurs inside the pores. The homogenisation technique recently proposed by Gao et al. (2016) was known to correctly predict the behaviour of standard acoustic foams. In this contribution, it has been used to obtain frequency dependent effective properties of the acoustic metafoam and to identify the parameters that have a dominant influence on the macroscopic behaviour of the material. It has been demonstrated that both the presence of the resonator and the fluid viscosity have a significant influence on the effective material properties and as a result on the sound attenuation behaviour. A good match between DNS solutions and homogenised results has been observed showing that the framework is able to capture the extraordinary behaviour of acoustic metafoams. However, it is worth mentioning that the extension of Biot's theory proposed by Gao et al. (2016) is required to capture the effect of local resonance. This work paves the way for designing realistic poro-elastic materials with enhanced performance at low frequency, since it enables modelling materials with complex microstructures, while taking into account their realistic sizes.  The macroscopic constitutive relations (14), (15), (16) and (17), using the Voigt notation, and taking into consideration the unit cell symmetries, are expressed as: