Computational design of locally resonant acoustic metamaterials

by   D. Roca, et al.
Universitat Politècnica de Catalunya

The so-called Locally Resonant Acoustic Metamaterials (LRAM) are considered for the design of specifically engineered devices capable of stopping waves from propagating in certain frequency regions (bandgaps), this making them applicable for acoustic insulation purposes. This fact has inspired the design of a new kind of lightweight acoustic insulation panels with the ability to attenuate noise sources in the low frequency range (below 5000 Hz) without requiring thick pieces of very dense materials. A design procedure based on different computational mechanics tools, namely, (1) a multiscale homogenization framework, (2) model order reduction strategies and (3) topological optimization procedures, is proposed. It aims at attenuating sound waves through the panel for a target set of resonance frequencies as well as maximizing the associated bandgaps. The resulting design's performance is later studied by introducing viscoelastic properties in the coating phase, in order to both analyse their effects on the overall design and account for more realistic behaviour. The study displays the emerging field of Computational Material Design (CMD) as a computational mechanics area with enormous potential for the design of metamaterial-based industrial acoustic parts.



There are no comments yet.


page 6

page 8

page 12

page 13


Structured model order reduction for vibro-acoustic problems using interpolation and balancing methods

Vibration and dissipation in vibro-acoustic systems can be assessed usin...

Acoustic Structure Inverse Design and Optimization Using Deep Learning

From ancient to modern times, acoustic structures have been used to cont...

Removal of spurious solutions encountered in Helmholtz scattering resonance computations in R^d

In this paper we consider a sorting scheme for the removal of spurious s...

Real-time 3-D Mapping with Estimating Acoustic Materials

This paper proposes a real-time system integrating an acoustic material ...

Material absorption-based carrier generation model for modeling optoelectronic devices

The generation rate of photocarriers in optoelectronic materials is comm...

Modeling the Repetition-based Recovering of Acoustic and Visual Sources with Dendritic Neurons

In natural auditory environments, acoustic signals originate from the te...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 counter-intuitive 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 low-frequency 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, LRAM-based 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 rubber-coated metal spheres embedded in polymer matrices, while more recently, Claeys et al. (2016) have carried out tests with a fully 3D-printed design with internal resonators capable of achieving also good attenuation properties in the low-frequency range. The idea of LRAM-based insulation panels has already been explored, both theoretically and experimentally, for membrane-type (Yang et al., 2010; Zhang et al., 2013) and plate-type (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 Bloch-Floquet 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 so-called 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 optimization-based 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 well-established 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 3D-printing 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 trial-and-error 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 low-frequency regime where the condition


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 Hill-Mandel 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:


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 Hill-Mandel principle, i.e.


along with kinematic restrictions that link the macroscopic displacements, , and strains, , with their microscale counterparts, and , namely


one can obtain a Lagrange-extended 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


An in-depth 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 Galerkin-based finite elements discretization, this so-called extended RVE system has the form:


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:

  • Quasi-static system ()

  • Inertial system ()


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:

  • The macroscopic strain accelerations are negligible, i.e. , allowing for the system (9) to actually behave quasi-statically, thus .

  • The RVE’s topology exhibits symmetry with respect to its geometric centre, which allows us to assume, for the subsystem (10), .

From the quasi-static part of the system, one can derive an expression for the macroscopic stress that reads



is an effective constitutive tensor.

On the other hand, from the inertial subsystem, it is possible to obtain the macroscopic inertial force as


where is the RVE average mass density and the second term in equation (12) represents the contribution of coupled micro-inertial 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)


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 stress-strain 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 rate-dependent effects, a more enriched Kelvin-Voigt model will be introduced (Krushynska et al. (2016); Lewinska et al. (2017)), so that the stress-strain relation becomes


where remains as the fourth-order constitutive tensor for an isotropic, linear elastic material and assumes the role of an analogous fourth-order viscous tensor. Since for most polymer-type materials (potential candidates as dissipative coating materials), rate-dependency affects mainly the deviatoric component of the strain velocity, typically the viscosity tensor will be considered as


where is the materials’ viscosity distribution and is the deviatoric fourth-order tensor, defined in index notation as


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


    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


    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 low-frequency 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 non-diagonal, 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 non-diagonal. However, the degree of coupling with non-relevant 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 level-set 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 lower-bound 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


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 level-set function such that


so that the function becomes the unknown of the problem in a variational context.

Figure 1: RVE configuration for the LRAM topology optimization (left) and typical effective LRAM density diagram depicting the frequency bandgap corresponding to negative densities (right). In order to ensure LRAM-like behaviour, a fixed matrix of material is considered, , with constant material properties and . Only in the inside region the characteristic function is allowed to change giving rise to dense material volumes corresponding to inclusions, , and soft material volumes corresponding to void/coating material, . The resulting resonance frequencies and associated to the restricted and unrestricted system modes, respectively, determine the lower and upper bounds of the target band, whose lower bound is matched to .

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


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 mass-normalized 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


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.


where in this case, the only relevant that correspond to upper bounds of the frequency bandgaps can be identified by


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.

Figure 2:

Global procedure to evaluate the transmission loss of a LRAM panel with a topology optimized design. The material distribution is obtained from the topology optimization algorithm. This result is used to build a RVE with the actual material properties from which the effective properties are computed by employing the multiscale homogenization framework. Finally, a macroscale analysis is performed over a slice of the panel, imposing displacement and traction compatibility conditions with the incoming and outgoing waves and periodic boundary conditions at the material boundaries (in order to simulate the infinite extension of the panel in the vertical direction). This analysis is performed in the frequency domain for several test frequencies allowing us to obtain the transmission and reflection coefficients of the panel,

and , respectively.

In this regard, the objective function proposed to minimize is given by




subject to the state-equations (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

s.t. (28)

A time-marching technique is used to update the problem with a pseudo-time variable, from an initial state (typically all the design area full with dense material) towards a problem’s solution. In this case, a Hamilton-Jacobi approach has been considered in which the problem’s evolution has been defined in a rate form as


allowing us to obtain the updated value of function from the previous iteration step through a straightforward time discretization of equation (29)


where is a pseudo-time 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 Hamilton-Jacobi 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 (non-design 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 non-relevant 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 horizontally-oriented 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
Table 1: Material properties (Calius et al., 2009).
Figure 3: Representation of the valid regions for achievable first target resonating frequencies in terms of the RVE size. Upper and lower limits and the so-called non-feasible region have been obtained with the material and numerical properties for this example, which are listed in Table 1. The green shaded area corresponds to the region of achievable first target resonating frequencies within the range of interest.

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 3-step 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 plane-strain 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

Figure 4: (a) Evolution of the objective function upon each iteration for the frequency fitting case (top) and the frequency fitting along with bandgap maximization (bottom). For , a trade-off must be met between the frequency fitting part and the bandgap maximization component, which makes it more difficult for the algorithm to converge. (b) Evolution of the volume fraction upon each iteration. Note that the colorbar represents the volume fraction distribution of each material in the domain. The fixed matrix (blue) represents 19% of the volume in both cases, and the remaining 81% is distributed between inclusion (grey) and coating phases (orange).
Figure 5: Evolution of the first relevant resonance frequencies for the restricted and unrestricted systems upon each iteration for the frequency fitting case (top) and frequency fitting along with bandgap maximization case (bottom). The topology determined by the level-set function at several iteration steps is also shown. Note the sudden jump in frequencies coincides with the iteration when the inclusion is disengaged. This gap gives an idea of the region of unattainable frequencies due to numerical issues.
Figure 6: First relevant resonance modes and frequencies for the restricted and unrestricted systems for the topologies resulting from the frequency fitting optimization (left) and the frequency fitting coupled with the bandgap maximization (right). The solid black lines represent the material interfaces of the non-deformed state.

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, non-relevant modes resulting from the modal analysis. A value of has been considered as a pseudo-timestep for the time-marching 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 full-material 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


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

Figure 7: Dispersion diagram for the topologies resulting from the frequency fitting optimization (left) and the frequency fitting coupled with the bandgap maximization (right). The shaded areas correspond to the bandgap regions corresponding to purely imaginary wavenumbers in the case without viscosity (). The units for the viscosity parameter are Pas. Note that the real part of the normalized wavenumber corresponds to the solid lines and ’x’ markers, while its imaginary component is represented through dashed lines and ’+’ markers in the same axes. The results obtained with the multiscale homogenization procedure are compared with those obtained applying Bloch-Floquet boundary conditions on the same RVEs, for validation purposes.

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 micro-inertial 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 non-null imaginary component of the wavenumber indicates the presence of a frequency bandgap. In Figure 7, the results of applying Bloch-Floquet 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 :


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


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


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

Figure 8: Transmission loss (TL) for the topologies resulting from the frequency fitting optimization (solid blue line) and the frequency fitting coupled with the bandgap maximization (red dashed line). The light-shaded arrows indicate the bandwidth of the effective attenuation bands corresponding to the first uninterrupted attenuations dB. The units for the viscosity parameter are Pas.

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 low-frequency 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 Bloch-Floquet 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 Bloch-Floquet 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 low-frequency range of interest where the local resonance of the RVE takes place.


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 DPI2017-85521-P 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


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


After a standard assembly process, denoted here with the big A symbol, one can obtain the global damping matrix


and the extended RVE dynamic system results


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,


where corresponds to the shape functions for the type of elements considered.

Since the hypotheses for the split of the system (41) into its quasi-static and inertial components still hold, despite the introduction of viscoelastic effects in the framework, one obtains:

  • Quasi-static system


    For the sake of generality, let us assume that kinematic conditions are imposed directly on the system through


    where , and is the centroid of the RVE (). In discretized form, this is




    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,


    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,


    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


    Taking the time derivative of equation (50) gives


    where the hypothesis for the quasi-static 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




    Finally, substituting expression (52) into equation (49), yields


    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 quasi-static 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


    Following a similar procedure than in the quasi-static 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


    which can be expressed in discretized form as




    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,


    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


    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)


    where here refer to the squared natural frequencies and are the associated mass-normalized vibration modes (i.e. and )111

    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 mass-normalized eigenmodes, and as the column vector with their corresponding modal amplitudes, one can express


    and equations (59) and (60) become


Appendix B Sensitivity of the LRAM topology optimization cost function

For any functional of the form


the VTD at any point in the domain is given by


where is the point-Dirac’s delta shifted to point fulfilling


On the other hand, the evolution of the cost function can be computed, accounting for equation (29), as


where equation (66) has been considered. From equation (20), it can be proven that


where stands for the line/surface-Dirac’s delta shifted to the zero level-set of in , i.e. , thus fulfilling


Replacing equations (69) and (70) into equation (68) yields


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


Now, in order to compute and let us take the derivative of the state-equations (21) and (23<