1 Motivation
The concept of Locally Resonant Acoustic Metamaterials (LRAM) has been object of growing interest among the scientific and technological communities in recent years. The notion of metamaterials emerged in the late 1990s as a new kind of engineered materials, by optimizing their morphology or arrangement at lower scales, capable of exhibiting properties “on demand” that are not found in naturally occurring materials. In the specific case of acoustic metamaterials, these properties involve exotic or counterintuitive behaviour caused by the interaction of the material with acoustic or elastic waves’ propagation features. These type of materials could be used, for instance, to design acoustic insulation panels that target specific frequency bands (especially in the lowfrequency range, i.e. below 5000 Hz, where most sources of environmental noise are produced). In contrast to conventional acoustic panels, which require either a large thickness or high mass density in order to provide effective sound attenuation at lower frequencies, LRAMbased panels can achieve good levels of transmission loss in the whole frequency range of interest with relatively thin, lightweight designs.
Scientific research in the field of metamaterials started in the late 19th century with the works of Floquet and Rayleigh among others, who studied phenomena related to the propagation of waves in periodic systems. However, it was not until the beginning of this century when the first implementation of an acoustic metamaterial capable of stopping waves from transmitting in a certain frequency band was reported by Liu et al. (2000). Later on, Ho et al. (2005) and Calius et al. (2009) performed other experimental demonstrations with silicone rubbercoated metal spheres embedded in polymer matrices, while more recently, Claeys et al. (2016) have carried out tests with a fully 3Dprinted design with internal resonators capable of achieving also good attenuation properties in the lowfrequency range. The idea of LRAMbased insulation panels has already been explored, both theoretically and experimentally, for membranetype (Yang et al., 2010; Zhang et al., 2013) and platetype (Badreddine Assouar et al., 2012) acoustic metamaterials.
On the other hand, computational models for the study and characterization of LRAMs have been traditionally focused on periodically repeated microstructures where BlochFloquet boundary conditions can be applied (Hussein, 2009; Hussein et al., 2014). More recent developments include the works of Sridhar et al. (2016, 2017) and Roca et al. (2018), in which computational homogenization frameworks accounting for inertial effects have been proposed, capable of capturing the local resonance phenomena that characterizes such kind of materials.
Aiming at trying to enlarge the effective attenuation band achieved by the local resonance phenomenon in metamaterials is a challenging task, since the socalled frequency bandgaps produced by typical LRAM designs, which are the source of their attenuation capabilities, tend to be too narrow in the frequency spectrum. In this regard, several proposals have been made recently in the literature in order to find a solution for this problem. For instance, Matsuki et al. (2014) proposed a topology optimizationbased method to obtain optimal LRAM configurations with multiple attenuation peaks, which can be considered one of the first attempts to apply optimization procedures to LRAM materials. Other approaches are focused on taking advantage of the viscoelastic nature of the coating materials in typical LRAM configurations. The first notions on the beneficial effects of viscous damping in acoustic metamaterials were reported by Hussein and Frazier (2013), who introduced the concept of metadamping to refer to the damping emergence phenomenon produced as a result of combining the effects of local resonance with viscous dissipation. The concept of metadamping has also been explored in more detail in subsequent works Frazier and Hussein (2015) and the idea of acoustic metamaterial configurations based on this (phononic resonators) has been proposed by DePauw et al. (2018). This phenomenon has also been studied more recently in the context of acoustic metamaterials in the work of Manimala and Sun (2014), where they show, by means of an analytical approach, the dispersion properties of LRAMs with different models of viscoelastic (dissipative) behaviour for the coating material, and the works of Krushynska et al. (2016) and Lewinska et al. (2017), in which generalized viscoelastic modelling is introduced in the study of the attenuation performance of LRAMs.
In this paper, the authors attempt at setting a computational based methodology for the modelling, analysis and design of metamaterial parts exhibiting local resonance phenomena by combining three wellestablished complementary computational tools: (1) a multiscale hierarchical homogenization procedure specifically devised for acoustic problems, described in Roca et al. (2018), (2) the exploitation of Reduced Order Modelling (ROM) techniques, to minimize the resulting computational cost, based on selective projections of the RVE behaviour onto the space spanned by its natural modes and (3) topology optimization techniques.
A design strategy is proposed to optimize the attenuation performance of an insulation panel made of LRAMs for a target band of frequencies, by optimizing the topology of the material at the mesoscale. This allows considering its industrial manufacture by means of emergent 3Dprinting or similar techniques. The proposed design strategy is based on (a) fitting the lower bound of the target band with some natural resonance frequency of the material at the design scale, and (b) maximizing the target band’s bandwidth. The resulting design exhibits acoustic insulation properties much improved, in comparison to those that could have been obtained by simple trialanderror procedures, which proves the benefits of the considered CMD methodology.
2 Multiscale modelling of LRAMs
The computational homogenization framework introduced in Roca et al. (2018) has been used here as the base model upon which the methodology for the design of LRAMs will be built. The framework can be applied to problems where a separation of scales is present, for instance, allowing us to identify a representative volume element (RVE), typically a unit cell, in a macroscopic structure. This separation of scales is established in terms of the macroscopic characteristic wavelength, , which has to be larger than the characteristic size of the lower scale, , otherwise the validity of the homogenization model cannot be guaranteed. This is not an impediment to deal with LRAMs considering they are designed to operate in the lowfrequency regime where the condition
(1) 
is easily satisfied. In fact, expression (1) is also a condition required for local resonance phenomena to arise (Hussein et al., 2014).
For clarity purposes, from now on magnitudes referring to the microscale will be subscripted by , in order to distinguish them from their macroscopic counterparts. Additionally, space coordinates for the macroscale will be referred by , while those for the microscale will be referred by , when necessary.
The framework is grounded on a generalization accounting for inertial effects of the classical HillMandel principle (Blanco et al., 2016), in which the macroscale is assumed to behave as a Cauchy’s continua, thus satisfying the classical postulates of linear and angular momentum:
(2)  
(3) 
where is the macroscopic stress, is the macroscopic inertial force, refers to the time variable, stands for the time derivative of , and refers to the spatial macroscopic domain. Note that body forces have not been considered, for the sake of simplicity, since they are not relevant in the context of acoustic problems that are tackled here.
By applying an energetic equivalence between scales, which is given by a variational statement of the generalized HillMandel principle, i.e.
(4) 
along with kinematic restrictions that link the macroscopic displacements, , and strains, , with their microscale counterparts, and , namely
(5) 
one can obtain a Lagrangeextended dynamic system of equations for the RVE in which the Lagrange multipliers, and respectively, corresponding to the reactions to the minimal kinematic restrictions given by equation (5) can be identified as
(6)  
(7) 
An indepth explanation of the theory and further details on the derivation of these terms can be found in Roca et al. (2018). Note that the angle brackets symbol is used to refer to the average volume integral over the RVE, i.e. .
After a Galerkinbased finite elements discretization, this socalled extended RVE system has the form:
(8) 
where and are the standard mass and stiffness matrices of the RVE system,
is the column vector with the nodal values for microscale displacement field,
, while and are, respectively, the Lagrange multipliers associated to the kinematic restrictions over displacements and their gradient, which are imposed by the matrices and , respectively.Note that, in the system (8), the displacement and strain of the associated point, , , become actions and, as indicated by equations (6) and (7), one can relate the resulting Lagrange multipliers, and , with the macroscopic inertial force and stress at that point, and , respectively. As it is shown in Roca et al. (2018), the system (8) can be split into:

Quasistatic system ()
(9) 
Inertial system ()
(10)
The split is based on considering each action separately in each subsystem along with certain hypotheses, which are suitable for the study of local resonance phenomena in acoustic problems. In particular:
From the quasistatic part of the system, one can derive an expression for the macroscopic stress that reads
(11) 
where
is an effective constitutive tensor.
On the other hand, from the inertial subsystem, it is possible to obtain the macroscopic inertial force as
(12) 
where is the RVE average mass density and the second term in equation (12) represents the contribution of coupled microinertial effects (through the matrix ) whose behaviour is dictated by the reduced set of uncoupled equations resulting from the modal projection of the inertial subsystem (10)
(13) 
where is a diagonal matrix containing the relevant natural frequencies of the inertial subsystem. More details on the derivation of these terms can be found in Roca et al. (2018).
3 Modelling the viscoelastic behaviour in LRAMs
The consideration of viscoelastic phenomena affects the model for the stressstrain relation. While in Roca et al. (2018), materials in the microscale where assumed to behave as linear elastic solids, here, in order to account for ratedependent effects, a more enriched KelvinVoigt model will be introduced (Krushynska et al. (2016); Lewinska et al. (2017)), so that the stressstrain relation becomes
(14) 
where remains as the fourthorder constitutive tensor for an isotropic, linear elastic material and assumes the role of an analogous fourthorder viscous tensor. Since for most polymertype materials (potential candidates as dissipative coating materials), ratedependency affects mainly the deviatoric component of the strain velocity, typically the viscosity tensor will be considered as
(15) 
where is the materials’ viscosity distribution and is the deviatoric fourthorder tensor, defined in index notation as
(16) 
with being Kronecker deltas (1 for and 0 otherwise).
Since the hypotheses for homogenization can be compatible with the introduction of this additional effect, the model still holds and the formulation naturally adapts to accommodate this new term. A detailed development of the model’s equations with this additional term is provided in A, so in this section, only the changes in the results will be discussed, i.e.:

The first relevant difference is in the expression for the macroscopic stress, which now reads
(17) Notice the appearance of an effective viscous tensor . This is not surprising considering the viscoelastic model assumed for the microscale in equation (14), which has an analogous form.

The second difference appears in the projected inertial system of reduced degrees of freedom. While formally the expression for the macroscopic inertia is the same than in equation (
12), an additional damping term appears in former equation (13), which now reads(18) The new matrix is responsible for damping the vibration near the resonance frequencies of the RVE. While this counteracts the effects of local resonance phenomena in the macroscale (especially the higher the frequency becomes), in a relatively lowfrequency regime (where LRAMs operate), and for certain levels of viscosity, it can provide the beneficial effect of extending the effective attenuation band. This phenomena will be observed and further described in section 5.
It is important to notice that, while the system of equations (13) is fully uncoupled, which allows us to perform an effective reduction of the number of degrees of freedom that need to be considered in the analysis, the same cannot be guaranteed for the system (18), due to the presence of the matrix , which is nondiagonal, in general, and can make the homogenization model more computationally expensive than in the case where no viscoelastic effects are considered. This increase in computational cost is related to the degree of copuling existing in the matrix , which ultimately depends on the RVE topology, material properties and modelling of viscoelastic effects in the microscale. For instance, the system (18) would remain fully uncoupled only in cases where the damping matrix is already diagonal or either it is proportional to the stiffness and/or the mass matrices (typically known as proportional Rayleigh damping model), which cause the matrix to be diagonal. In the specific cases accounting for viscous effects considered in this work, strain rate dependence is only considered in one of the material components, i.e. the coating, thus, in general, the matrix is nondiagonal. However, the degree of coupling with nonrelevant modes will be small enough to be negligible in practice..
4 Topological design of LRAMs
Aiming at obtaining LRAMs with topologies designed to achieve better attenuation properties, i.e. increasing the levels and range of transmission loss (in the specific case of acoustic insulation panels), a levelset based topology optimization strategy is proposed here with two main objectives:

Fit the relevant resonance frequencies of the resulting topology into a targeted band. This is done by matching the lowerbound of the target band with a relevant resonance frequency of the RVE (see Figure 1).

Maximize the bandwidth of the target band in terms of the topology of the RVE materials.
In such methodologies, a suitable cost function is minimized, in a variational way, with respect to a characteristic function
(19) 
that defines the material distribution in the design domain , taking values of 1 for dense material regions (inclusions): and 0 for soft material regions (coating/void): (). These regions are typically determined by a smooth levelset function such that
(20) 
so that the function becomes the unknown of the problem in a variational context.
Prior to presenting the optimization problem itself, let us first review the RVE equations that will be required to solve this problem. First, in order to obtain the RVE’s relevant resonance frequencies, the modal problem of the restricted system must be solved
(21) 
where and are the resulting stiffness and mass matrices once the kinematic restrictions on the microfluctuation field have been applied. In this specific case, since the local resonance phenomenon occurs at frequencies corresponding to internal vibration modes, a good approximation to meet our goals consists of prescribing all RVE boundaries. The terms and correspond to the squared natural frequencies and massnormalized vibration modes of the system. According to equation (12), only those modes such that their corresponding columns in the generalized coupling matrix are larger than a certain tolerance , should be considered as relevant. Thus, the set of relevant resonance frequencies is given by
(22) 
Furthermore, since transmission loss peaks in LRAM panels are closely related to frequency bandgaps, whose lower and upper limits can be identified, respectively, with the relevant resonance frequencies of the restricted and unrestricted RVE system, and , as reported in Roca et al. (2018) (see Figure 1 for a schematic representation), one can maximize their bandwidth, for instance, by minimizing the ratio , where comes from the modal problem considering the mass and stiffness matrices of the RVE system prior to applying the kinematic restrictions, i.e.
(23) 
where in this case, the only relevant that correspond to upper bounds of the frequency bandgaps can be identified by
(24) 
Note that, in this case, the additional condition needs to be applied in order to avoid rigid body translation modes, which are relevant according to condition , but do not define the upper bounds of any bandgap.
In this regard, the objective function proposed to minimize is given by
(25) 
with
(26)  
(27) 
subject to the stateequations (21) and (23). In equations (25) to (27), is a weighting parameter to establish the relative importance of each term in the global objective function, is the imposed targeted squared frequency to fit with the first relevant squared resonance frequency for the restricted RVE system , while refers to the first relevant squared resonance frequency for the unrestricted RVE system. Note that, as long as , the objective function will be bounded , so that it reaches its minimum value when all the desired objectives are fulfilled. For the sake of simplicity, only the first frequency band is targeted to fit in this framework, but equation (25) can be easily extended to target multiple frequency bands simply by adding the corresponding terms in the cost function.
The optimization problem then reads
FIND:  
FULFILLING:  
s.t.  (28)  
A timemarching technique is used to update the problem with a pseudotime variable, from an initial state (typically all the design area full with dense material) towards a problem’s solution. In this case, a HamiltonJacobi approach has been considered in which the problem’s evolution has been defined in a rate form as
(29) 
allowing us to obtain the updated value of function from the previous iteration step through a straightforward time discretization of equation (29)
(30) 
where is a pseudotime step, is a parameter and is the topological sensitivity, evaluated at point , of the cost function (25), here named as the Variational Topological Derivative (VTD) of the functional . In B it is proven that the iterative scheme in equation (30) yields an iterative descend of the cost function , i.e. , a crucial aspect for the convergence of the HamiltonJacobi algorithm.
To ensure that local resonance phenomena arise in the computed designs throughout the optimization process, the frequency fitting will always be required, and typically aimed at the smallest relevant resonance frequency achievable (i.e. ) which, for a given set of material properties, will be constrained by the dimensions of the design domain (in this case the RVE). Furthermore, the RVE will consist of a fixed matrix material frame (nondesign domain) so that the actual design domain is the inclusion/coating distribution on the inside (see Figure 1).
Some additional hypothesis have been made in order to avoid spurious modes resulting from the modal analysis, which helps the optimization algorithm to become more stable and converge to the desired solutions. These hypothesis are listed and explained below:

The matrix fixed frame is considered infinitely stiff so that no deformation modes, which are nonrelevant in this context, appear in the modal analysis. This is done to prevent modes in which the matrix interacts with the other materials, especially in early steps of the algorithm, when the whole domain is filled with material.

The void/coating material is considered massless. Since the modes causing the local resonance phenomena to arise are those in which the inclusions vibrate inside the coating phase, the relevant properties that are to be considered will be the density of the inclusion phase and the stiffness of the coating material. Forcing the density of the coating material to zero (or a very small tolerance value), avoids the appearance of spurious modes in the modal analysis which greatly helps the identification of the relevant resonance modes.

Since the focus in this context is in horizontallyoriented modes (the panel will be subjected to plane waves propagating on the horizontal direction), all vertical degrees of freedom are prescribed in both the restricted and unrestricted systems. By restricting the analysis to a single dimension, both the identification of the relevant modes and the pairing of each bandgap limiting frequencies become much easier.
Material  Density  Bulk mod.  Shear mod. 

(kg/m)  (MPa)  (MPa)  
Epoxy  1180  
Steel  7780  
Silicone rubber  1300  0.63  0.04 
5 Application to the design of an acoustic insulation panel
Let us consider an infinitely large flat panel with a given thickness built with stacked LRAM unit cells (of size ) consisting of 3 material phases: an epoxy matrix frame at the boundaries along with a certain distribution of steel inclusions embedded in a silicone rubber coating (see Figure 1 for a graphical depiction of a typical RVE configuration). The material properties used in the examples that follow are listed in Table 1. The coating material will be considered viscoelastic, with the viscosity left out as a parameter in order to enable the possibility of evaluating the LRAM behaviour for various degrees of dissipation. Since the aim of this analysis is to assess the attenuation of acoustic waves through a slab of a designed LRAM panel, the transmission loss in the specific frequency range of interest (in this case below 3000 Hz) will be computed. To do so, a 3step analysis is performed. First, the topology optimization procedure explained in section 4 is used to obtain a material distribution for the LRAM design that meets the desired properties. Then, a RVE is built upon the results obtained from the previous step with the actual material properties, so the homogenization procedure detailed in sections 2 and 3 can be applied to compute the effective material properties associated to that LRAM design. Finally, an analysis on the macroscale is performed by computing the transmission loss in the desired frequency range for a flat panel composed of the homogenized LRAM subjected to acoustic plane waves. Details on each of these procedures is given in the following sections, while a summary of the global scheme is depicted in Figure 2. A planestrain 2D approach is considered in all the examples instead of a 3D setting simply to avoid the unnecessary complexity associated with them which, at least regarding the effects and conclusions that are expected to point out in this work, do not give any relevant additional insights. Therefore, for the sake of clarity in terms of interpretation of the results, the examples shown here are all 2D, even though both the formulation and the conclusions that can be extracted can be extended to 3D.
5.1 Topology optimization of the LRAM
A 2D structured mesh of quadrilateral elements with 4 Gaussian integration points has been used for the computations in this stage. The elastic properties of the matrix phase have been scaled by a factor to guarantee its behaviour as a rigid component, while the density of the coating has been scaled by to avoid spurious, nonrelevant modes resulting from the modal analysis. A value of has been considered as a pseudotimestep for the timemarching algorithm, along with an initial that makes all the design domain to be full of inclusion dense material (see equation (20)). The size of the RVE has been chosen to be 1
1 cm to ensure that the resulting resonating frequencies lie on the desired range. This is important since for a given set of material properties and a domain size, there are limitations in the achievable resonant frequencies. An obvious first upper limit can be found in the first resonant frequency obtained by the initial fullmaterial configuration, which depends on the properties of the inclusion phase as well as the RVE size and is typically well above (by several orders of magnitude) the desired frequency range. On the other hand, a theoretical lower limit can also be estimated as
(31) 
where the index here refers to each material. Furthermore, given the high contrast between the properties of the inclusion and coating components, there is a range of frequencies between these lower and upper limits for which the algorithm finds it more difficult to converge to a solution. This is because solutions in this range typically contain physically unstable solutions characterised by large changes in the resonance frequencies for small perturbations in the topology (in this case, caused by the appearance of unrealistically thin strings of material). Therefore, designs inside this zone should be avoided to preclude such unstable behaviour. Note also that these solutions can be easily avoided for a given target frequency by reducing the size of the design domain (see Figure 3 as an example).
For this example, a target frequency of Hz is selected and two different weighting parameters have been tested: one that causes only the frequency fitting to be considered in the objective function (), and another one where the weight is distributed equally among the fitting and the bandgap maximization parts (). Figure 4 shows the evolution of the objective function and the volume fraction of material, respectively. The first relevant resonating frequencies for the restricted and unrestricted problems can be seen in Figure 5.
These results show how the algorithm removes inclusion material in the topology, reducing the value returned by the objective function, as seen Figure 4 (a), and thus approaching the desired targets with each iteration step. It is worth noting the sudden jump produced around the 100th iteration in both cases. It is caused by the disengagement of the inclusion material from itself, which makes the coating phase to fully envelop it. Since the stiffness of the coating material is several orders of magnitude lower, the resulting configuration makes it easier for the enclosed inclusion to vibrate, which translates into a much lower resonance frequency.
Interestingly, while the volume fraction of material ends up being very similar in both cases, as seen in Figure 4 (b), the resulting topologies are quite different. Note also that the resonance frequency for the restricted problem follows a similar evolution in both cases (eventually meeting the target frequency of 1000 Hz), but the differences between their associated topologies make their respective unrestricted system’s resonance frequency to evolve differently (see Figure 5). In particular, when the bandgap maximization is also part of the objective function the algorithm tends to remove material from the matrix internal borders, concentrating the maximum amount of mass onto the central inclusion. This translates into a significant increase of the bandgap size, for a similar volume fraction of material distribution, by almost 6 times (around 600 Hz in the first case and up to 3500 Hz in the second).
It must be considered that the frequencies obtained from the optimization process will not be equal to those computed when actual material properties (without scaling factors) are considered, as it will be seen in the results of the next section. However, this is not a major issue for our purposes, since with the actual material properties, the system tends to be less restricted and so the resulting frequencies are smaller, which is often desirable. In any case, the size of the RVE can be adjusted accordingly to match the target frequencies, if desired.
5.2 Homogenization of the optimized LRAM
With the material distributions obtained in the previous step, RVEs are built, maintaining their size but with the actual material properties. The homogenization framework detailed in sections 2 and 3 is applied to compute the effective material parameters associated to each LRAM design. These include the effective constitutive tensor, , the average density , the coupling matrices related to microinertial phenomena, , , and the additional terms , , accounting for viscous effects. For these computations, 2D meshes of triangle elements with an average relative size of 0.02 and 3 Gaussian integration points have been used, which provide around 7000 degrees of freedom. Figure 6 shows the resulting first resonance modes and frequencies which, as previously anticipated in Remark 5.1, are smaller than the computed as part of the optimization process in the previous section. Note, however, that the predicted bandgap in the case with is much larger than in the other case (around 1500 to 450 Hz, which is between 3 and 4 times bigger). This can also be seen in the dispersion diagrams of Figure 7, which are calculated assuming a plane wave travelling in the direction in an infinite extension of the homogenized material domain (see Roca et al. (2018) for details on this computation). For the case without viscosity (), a nonnull imaginary component of the wavenumber indicates the presence of a frequency bandgap. In Figure 7, the results of applying BlochFloquet boundary conditions on the same RVEs (Hussein et al., 2014) are also shown to validate the proposed model. Note that there is good agreement between both computations, even for the cases where viscosity is present.
5.3 Transmission loss computation
Finally a macroscale analysis is performed over a homogenized panel from which the transmission coefficient is obtained. For the sake of simplicity, the panel is assumed to be surrounded by air (with density and sound propagation speed ) where a plane wave travels perpendicular to the panel’s surface. In order to find the transmission coefficient in this case, compatibility conditions for the normal component of the displacement and traction forces are imposed at the panel’s interfaces with the air media. In this context, one can assume the air at the left hand side to propagate as a wave at a certain frequency :
(32)  
(33) 
where and are the horizontal displacement and pressure fields, respectively, is the wavenumber and is the reflection coefficient, i.e. the fraction of the incident wave’s amplitude that is reflected at the panel’s interface. At the right hand side, only part of the wave is transmitted, so
(34)  
(35) 
where is the transmission coefficient, i.e. the fraction of the incident wave’s amplitude that is transmitted through the panel. To find and , a simple 2D FE discretization of quadrilateral elements with 4 gaussian integration points over a portion of the panel is performed. The resulting system is solved in the frequency domain for each frequency in the desired range, in this case from 0 to 3000 Hz. Periodic boundary conditions are applied on the top and bottom boundary nodes (to simulate the infinite extension of the panel in the vertical direction) and the compatibility conditions given by equations (32) to (35) are applied on the left and right boundary nodes. The resulting system can be expressed in matrix form as
(36) 
Since both and are functions of the reflection and transmission coefficients (as a result of applying the compatibility conditions), equation (36) can be reduced to a system with only and as unknowns (details on the derivation of this system are given in C). Thus, for a given frequency and viscosity , once the transmission coefficient is solved, the transmission loss is computed using the expression
(37) 
The results are shown in Figure 8 where the transmission loss has been computed for a panel with thickness of 1 cm (equivalent to a single RVE) and for both the topologies without and with the bandgap maximized. One can observe that the effect of the bandgap maximization () translates also into a larger frequency range of effective attenuation. For instance, for attenuation levels above 40 dB, the effective bands in both cases start around 250 Hz, but they extend to 1180 Hz in the first case and to 1840 in the second (660 Hz increase).
It can also be observed that the viscosity has the beneficial effect of smoothing the undesired inverted resonance peaks but at the same time it also dampens the attenuation peaks caused by the local resonance phenomena occurring at the microscale. Depending on the desired attenuation performance, one can take advantage of the viscoelastic properties of the coating component to slightly increase the effective attenuation band or even bypass the undesired inverted resonance peak while keeping a continuous effective attenuation band with a good level of attenuation in the low frequency range of interest (as it can be seen for the case of Pas, where viscous effects are more relevant).
6 Concluding remarks
While the methodology presented has been used in the study of flat panels under plane acoustic pressure waves, it can be suitably adapted to more complex cases still satisfying the lowfrequency range restriction. On the one hand, the objective function of the topology optimization algorithm can be adapted to tackle one or multiple frequencies depending on the user interest. On the other hand, the ability of the homogenization procedure to provide a set of effective constitutive properties makes it possible to easily study the behaviour of complex geometries in the macroscale under several sets of boundary conditions and external actions, including, for instance, the characterization of transient states or multiple layer configurations, so long as they satisfy the hypotheses considered. This is in contrast to other homogenization approaches, such as those based on the BlochFloquet theory, which have the ability to study the effective behaviour of RVEs without restrictions in terms of frequencies, but that are limited to predetermined macroscopic configurations, namely, infinitely periodic structures, in the case of the BlochFloquet theory, that cannot account for the introduction of other more realistic sets of boundary conditions.
In this regard, the design procedure presented in this work offers a powerful tool for the design of LRAM acoustic insulation devices tackling the low frequency range. By combining the optimization algorithm for designing the RVE topology with the homogenization method for characterising the material’s performance, one can easily assess the transmission loss of a panel for a broad range of frequencies. The studies carried out show the influence of the RVE topology in achieving certain properties in the final design such as, for instance, an increased bandgap size for a given frequency range with the same amount of material, and thus the same mass. It is worth noting that while the target frequencies and choice of materials greatly determine the size of the RVE, its topology optimization can be used to reduce the effective density making the resulting panel suitable for lightweight applications. On the other hand, it has been observed that the presence of highly viscoelastic materials in the design generally affects negatively the local resonance performance of the RVE. However, in certain viscosity ranges, one can take advantage of the resulting damping effects to smooth undesired resonance peaks in the macroscale and even join attenuation bands, specially in higher frequency ranges, while still maintaining good attenuation levels in the lowfrequency range of interest where the local resonance of the RVE takes place.
Acknowledgements
This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Proof of Concept Grant agreement nº 779611) through the project “Computational catalog of multiscale materials: a plugin library for industrial finite element codes” (CATALOG). The authors also acknowledge the funding received from the Spanish Ministry of Economy and Competitiveness through the research grant DPI201785521P for the project “Computational design of Acoustic and Mechanical Metamaterials” (METAMAT).
Appendix A Introduction of viscoelastic effects in the homogenization framework
Considering equation (14) for the stress definition in the microscale, causes a damping matrix, , to arise from the Finite Element discretization of the RVE system, in addition to the standard mass and stiffness matrices, and . This damping matrix appears only in the elements where viscoelastic effects are considered, in which it is defined as
(38) 
where gives the derivatives of the shape functions for the type of elements considered, and is the element’s viscous tensor, which can be expressed, following Voigt’s notation, in terms of the associated material’s viscosity, , as
(39) 
After a standard assembly process, denoted here with the big A symbol, one can obtain the global damping matrix
(40) 
and the extended RVE dynamic system results
(41) 
Note also that matrices denoted by and appear in system (41). These come from the discretization of the minimal kinematic restrictions, as explained in detail in Roca et al. (2018). As a reminder,
(42)  
(43) 
where corresponds to the shape functions for the type of elements considered.
Since the hypotheses for the split of the system (41) into its quasistatic and inertial components still hold, despite the introduction of viscoelastic effects in the framework, one obtains:

Quasistatic system
(44) For the sake of generality, let us assume that kinematic conditions are imposed directly on the system through
(45) where , and is the centroid of the RVE (). In discretized form, this is
(46) with
(47) and being a matrix that imposes the desired boundary conditions on the discretized microfluctuation field, . Typically, for the kind of problems tackled, periodic boundary conditions along with prescription of certain degrees of freedom to prevent rigid body motions offer good results. In such cases,
(48) where superscripts 0, , and are used to refer to prescribed, internal, and periodic boundaries degrees of freedom, respectively.
Note that, from the definition of matrix (see equation (43)), it can be seen that . Thus, by premultiplying the first equation in system (44) by and introducing expression (46) allows us to find, after some algebraic manipulation,
(49) Assuming the kinematic restriction imposed is compatible with that associated to the Lagrange multiplier, projecting the system into a set satisfying the microfluctuation field boundary conditions, i.e. premultiplying the first equation in system (44) by , results in
(50) Taking the time derivative of equation (50) gives
(51) where the hypothesis for the quasistatic system, (which makes ), has been considered.
The solution for obtained from equation (51) in terms of the macroscopic strain rate, , can be introduced into equation (50) to compute the solution field
(52) with
(53) Finally, substituting expression (52) into equation (49), yields
(54) Note that in equation (54) assumes the role of an effective constitutive tensor and it is the same than the one obtained in Roca et al. (2018), without considering viscoelastic effects. In fact, the quasistatic influence of viscoelasticity in the model is accounted by the new term , which acts as a sort of effective viscous tensor relating the macroscopic stress with strain rates.

Inertial system
(55) Following a similar procedure than in the quasistatic case, in order to maintain the generality of the formulation, the specific kinematic conditions considered here will also be introduced directly in the system (55) by taking
(56) which can be expressed in discretized form as
(57) with
(58) and being, again, a matrix imposing the desired boundary conditions over the microfluctuation field (see, for instance, eq. (48)).
In this case, it is easy to verify, from the matrix definition in equation (42), that . Also note that, since is a rigid body translation mode, it belongs to both the kernels of and (i.e. and ). These properties allow us to obtain, after premultiplying the first equation in system (55) by and some algebraic manipulation,
(59) Again, projecting the system (55) into a set satisfying the imposed kinematic restrictions on the microfluctuation field through (nullifying the effect of the Lagrange multiplier), allows us to obtain
(60) In order to make the model computationally efficient, a model order reduction is performed by projecting the solution field onto the space spanned by the eigenmodes of the undamped system (60)
(61) where here refer to the squared natural frequencies and are the associated massnormalized vibration modes (i.e. and )^{1}^{1}1
Code functionalities allowing the computation of only the smallest eigenvalues in magnitude and their associated eigenvectors (which are potential candidates to become relevant modes in the frequency range of interest) have been used, this translating into a reduction of the computational cost of the modal problem evaluation.
. Now, defining as the diagonal matrix with only the natural frequencies in the range of interest, as a matrix that contains their associated massnormalized eigenmodes, and as the column vector with their corresponding modal amplitudes, one can express(62) and equations (59) and (60) become
(63) (64)
Appendix B Sensitivity of the LRAM topology optimization cost function
For any functional of the form
(65) 
the VTD at any point in the domain is given by
(66) 
where is the pointDirac’s delta shifted to point fulfilling
(67) 
On the other hand, the evolution of the cost function can be computed, accounting for equation (29), as
(68) 
where equation (66) has been considered. From equation (20), it can be proven that
(69) 
where stands for the line/surfaceDirac’s delta shifted to the zero levelset of in , i.e. , thus fulfilling
(70) 
Replacing equations (69) and (70) into equation (68) yields
(71) 
Equation (71) proofs the descending character of the cost function along time/iteration evolution.
The sensitivity of the cost function (25
) will be evaluated using the variational topological derivative (VTD). To do so, first the chain rule will be applied to equation (
25) to obtain(72) 
Now, in order to compute and let us take the derivative of the stateequations (21) and (23<
Comments
There are no comments yet.