Analysis of Low Frequency Acoustic Stop Bands in Cubic Arrays of Thick Spherical Shells With Holes

Dupont, Guillaume; Movchan, Alexander; Enoch, Stefan; Guenneau, Sébastien

doi:10.3389/fmats.2019.00050

ORIGINAL RESEARCH article

Front. Mater., 11 April 2019
Sec. Mechanics of Materials
Volume 6 - 2019 | https://doi.org/10.3389/fmats.2019.00050

Analysis of Low Frequency Acoustic Stop Bands in Cubic Arrays of Thick Spherical Shells With Holes

Guillaume Dupont¹

Alexander Movchan²

Stefan Enoch³

Sébastien Guenneau³^*

¹Aix Marseille Univ, CNRS, Centrale Marseille, IRPHE, Marseille, France
²Department of Mathematical Sciences, The University of Liverpool, Liverpool, United Kingdom
³Aix Marseille Univ, CNRS, Centrale Marseille, Institut Fresnel, Marseille, France

We analyse the propagation of airborne pressure waves through a three-dimensional array of rigid coated spheres (shells) in air. When we dig a channel terminated by an air cavity in each rigid shell we observe the appearance of a low frequency stop band. Each shell with a hole acts as a Helmholtz resonator supporting a low frequency localized mode. Isofrequency surfaces and contours reveal the strong anisotropy of the periodic structure at the edge of the stop band. A simple mechanical model of springs and masses allows for asymptotic estimates of the low frequency stop band for elongated channels. Increasing the radius of an air channel shifts up the position, and enlarges, the low frequency stop band. Adding holes in shells also shifts up the frequency of the stop band, and embedded shells lead to additional stop bands. Localization effect induced by a large defect in a periodic macrocell of Helmholtz resonators is finally investigated.

1. Introduction: Acoustic Metamaterials

In the tracks of photonic crystals, phononic crystals (Dowling, 2008) have provided a fillip for research in acoustic stop band structures (Kushwaha et al., 1993; Kafesaki and Economou, 1999) within which light or sound is prohibited to propagate due to multiple scattering between periodically spaced inclusions. In 2000, Liu et al. provided the first numerical and experimental evidence of frequency dispersive elastic parameters of locally resonant structures for elastic waves in three-dimensional arrays of thin coated spheres (Liu et al., 2000): The effective parameters were shown to turn negative where low frequency stop bands occur. This important work paved the way toward acoustic analogs of electromagnetic meta-materials, such as fluid-solid composites for the control of pressure waves (Auriault and Bonnet, 1985; Auriault, 1994; Liu et al., 2000; Goffaux et al., 2002; Ho et al., 2003; Elford et al., 2004; Hirsekorn et al., 2004; Li and Chan, 2004; Movchan and Guenneau, 2004; Fang et al., 2006; Mei et al., 2006; Guenneau et al., 2007b; Hou et al., 2007; Wang et al., 2008; Chalmers et al., 2009; Norris, 2009; Lemoult et al., 2011; Auriault and Boutin, 2012; Boutin, 2013; Boutin and Becot, 2015). Using asymptotic methods for fields in multi-structures (Kozlov et al., 1999; Movchan et al., 2002), it has been proposed to use arrays of cylinders with a split ring cross section as building blocks for two-dimensional localized resonant acoustic structures displaying negative refraction (Movchan and Guenneau, 2004; Guenneau et al., 2007a), based on analogies with split ring resonators (SRRs), introduced by Pendry in the context of electromagnetic waves almost 20 years ago (Pendry et al., 1999). Magnetic activity of metamaterials occurs near resonances of SRRs, which was essential for instance in the design of the first electromagnetic cloak (Schurig et al., 2006), and similar designs have been proposed for acoustic cloaks (Craster et al., 2013).

In a series of articles, the research group of Auriault developed asymptotic models of locally resonant structures (Auriault and Bonnet, 1985; Auriault, 1994; Auriault and Boutin, 2012), some of which predate the birth of acoustic metamaterials with the seminal work of the research group of Ping Shen (Liu et al., 2000). Low frequency stop bands in arrays of 2D Helmholtz resonators with elongated necks have been studied in (Movchan and Guenneau, 2004; Guenneau et al., 2007a) as mentioned above, with further predictions and measurements of sound transmission in elastic shells in air in Krynkin et al. (2010). However, 3D Helmholtz resonators with elongated neck position inside the cavity and periodically distributed in a fluid have been introduced in (Boutin, 2013; Boutin and Becot, 2015). Therein, Boutin and coauthors have shown that the macroscopic fields (pressure P and velocity V) satisfy the following equation of mass conservation and generalized Darcy Law:

\begin{array}{l} \nabla \cdot V = \frac{i ω Φ}{K} P + \frac{Q}{∣ Y ∣}, V = \frac{- T}{η} \nabla P & (1) \end{array}

where K is the effective bulk modulus (including thermal effects) of the matrix, ω the resonant frequency, i² = −1, Φ and T are the porosity and dynamic permeability tensor (including viscosity effects) when considering the resonators as perfectly rigid spheres, η is the viscosity of the fluid-matrix, ∣Y∣ is the volume of the 3D unit cell Y and Q is the flux emitted by the resonator into the fluid matrix in response to the pressure acting on it. The case of the Helmholtz resonator with a single elongated neck has been not only derived, but also experimentally tested, in Boutin (2013) and Boutin and Becot (2015). And the fluid need not be inviscid.

There is a vast amount of literature on the role of the neck in such Helmholtz resonators (Groeneweg, 1969; Alster, 1972; Orris et al., 1974; Gaunaurd and Uberall, 1982; Lim et al., 1990; Hinders et al., 1995; Selamet and Dickey, 1995; Mead, 1996; Aberg et al., 1997; Selamet et al., 1997; Chen et al., 1998; Baird et al., 1999; Seo et al., 2005; Ivansson, 2006; Duan et al., 2007; Zhao et al., 2007; Zhou et al., 2010; Wang and Mak, 2012; Li et al., 2013). Recently, periodic structures consisting of SRRs have been revisited thanks to refined homogenization techniques, which allow to rigorously link the geometrical parameters to frequency dependent effective parameters (Mercier et al., 2017). Analogs of SRRs have been proposed for in-plane elastic waves (Guenneau et al., 2007b), that allows for bending and rotational vibrations having lower resonant frequencies than compressional vibrations, and thus make good candidates for building blocks of mechanical metamaterials, which might make possible some mechanical cloak, see Kadic et al. (2013) for a review of metamaterials beyond electromagnetics. Milton, Briane and Willis provided a thorough mathematical frame for cloaking for certain types of elastodynamic waves in structural mechanics, in the framework of modified Willis equations (Milton et al., 2006). On the other hand, coupled in-plane pressure and shear waves were numerically shown to be detoured around a finite size obstacle by a specially designed cloak with an anisotropic heterogeneous elasticity tensor without the minor symmetries (Brun et al., 2009) without resorting to Willis media. Nonetheless, SRRs like in Guenneau et al. (2007b) make an interesting candidate to achieve such Cosserat-type metamaterials. Actually, a similar type of resonant elastic structure already led to a negatively refracting medium (Zhu et al., 2014), where the elastic chirality (which is a hallmark of a Cosserat medium) was put forward in conclusive experiments. The need for mechanical metamaterials with elasticity tensor without the minor symmetry motivates the extension of homgenization results of Mercier et al. (2017) to the Navier system. Importantly, the acoustic wave equation is invariant under coordinate changes (Norris, 2008, 2009), so acoustic metamaterials via geometric transform can thus achieve enhanced control of pressure waves via a simpler route. A simple example of a periodic system is a periodic lattice, whose dynamic response is well defined by Green's functions. Dynamic Green's functions in periodic lattices and their asymptotics were analyzed in Movchan and Slepyan (2014) and Vanel et al. (2016). Localization of electromagnetic waves within a two-dimensional grating of spheres embedded in the three-dimensional space and the high-frequency homogenisation approximations were studied in Maling et al. (2017).

A downfall of SRRs is that they do not allow for doubly negative acoustic parameters, so one would need a second type of inclusion within the periodic cell to achieve a goal of a negative refractive index in order to design a super lens like Pendry proposed in his seminal paper (Pendry, 2000). In 2004, Li and Chan proposed a design of negative acoustic metamaterial based on a multiple scattering theory approach (Li and Chan, 2004), which has been revisited recently by applied mathematicians using homogenization theory for bubbly fluids (Ammari et al., 2017a,b) with hybridization of Minnaert resonances (Minnaert, 1933). This is somewhat related to homogenization of high-contrast periodic structures (Figotin and Kuchment, 1998). Potential applications of bubbly media in super lensing and total absorption have been experimentally shown (Lanoy et al., 2015; Leroy et al., 2015). Such soft metamaterials (Brunet et al., 2015) are currently investigated by many groupings in the world, but they might not always be compatible with industrial processes.

Although SRRs might not achieve double negative acoustic parameters, they seem to be quite straightforward to engineer for ultrasonic waves. Fang et al. experimentally demonstrated a dynamic effective negative stiffness in a chain of air filled Helmholtz's resonators for ultrasonic waves (Fang et al., 2006). Moreover, it has been also shown using asymptotic techniques that surface water waves propagating within an array of fluid filled SRRs display a negative effective density (Farhat et al., 2009). A focussing effect through a finite array of such resonators was numerically achieved, with a resolution of a third of the wavelength. The negative effective density, reinterpreted as a negative effective gravity, has been further confirmed theoretically and experimentally (Hu et al., 2011, 2017). Fluid filled SRRs have been experimentally tested in a 17 meter long water channel for filtering effects with a dike using grooved vertical cylinders (Dupont et al., 2017). Similar filtering effects have been experimentally demonstrated for sound waves interacting with a doubly periodic array of Helmholtz's resonators shaped as soda cans (Lemoult et al., 2011). However, we note that in that case, there is no neck, just a hole, as soda cans are thin shells. In the present paper, we would like to revisit the concepts of SRRs in the case of pressure waves propagating in a three-dimensional array of 3D Helmholtz resonators with elongated necks. In order to simplify the mathematical setup, we shall consider airborne acoustic waves, in which case rigid (Neumann type) boundary conditions can be considered on resonant elements. Compared to problems of linear elasticity, the present study does not deal with dynamic degeneracies at low frequencies, which may occur for certain type of geometries of elastic systems, resulting in a group of very small eigenvalues being separated from the remaining spectrum.

2. Motivation: Spectral properties of a periodic array of rigid spheres

Let us start with an illustrative numerical result for a spectral problem for the Helmholtz operator within a periodic cubic array of rigid spheres: the unknown is a pressure wave field, here sound in air (wave speed c = 340m.s⁻¹). Neumann boundary conditions are prescribed on the contour of each defect and standard Floquet-Bloch conditions are set on an elementary cell of the periodic structure. The finite element formulation was implemented in the COMSOL Multiphysics Package to compute the eigenvalues and to generate the corresponding eigenfields. We present in Figure 1 the periodic structure we want to study and in Figure 2 the corresponding dispersion diagram for eigenfrequencies ω (in unit of rad.s⁻¹) as a function of the Floquet-Bloch parameter k (in unit of m⁻¹): along the horizontal axis we have the values of modulus of k, where k stands for the position vector of a point on the contour ΓXMU within the irreducible Brillouin zone. We note the absence of bandgaps with the presence of rigid spheres. This lack of intervals of forbidden frequencies motivates the present study: how can one create a stop band without further increasing the size of the rigid spheres?

FIGURE 1

Figure 1. Three-dimensionnal phononic crystal, periodic cell Y with a rigid sphere Ω in physical space and irreducible Brillouin zone ΓXMU of the periodic cell Y* in reciprocal space with the three components of the Floquet-Bloch vector k = (k_x, k_y, k_z).

FIGURE 2

Figure 2. Dispersion diagram for a periodic array (pitch d = 1 m) of spherical rigid inclusions (R = 0.4 m) representing the frequency ω (rad.s⁻¹) of pressure waves in air, vs. the wavenumber |k| (m⁻¹), projection of the Bloch vector k along the edges of ΓXMU.

3. Setup of the spectral problem: The continuum model

Let us first recall the finite element set-up. We consider the Helmholtz equation:

\begin{array}{l} \nabla . (\frac{1}{ρ (x, y, z)} \nabla p (x, y, z)) + \frac{ω^{2}}{K (x, y, z)} p (x, y, z) = 0 & (2) \end{array}

where ρ (kg.m⁻³), K (Pa) are the density and bulk modulus of the medium and ω (rad.s⁻¹) is the angular wave frequency of the pressure field p.

Due to the periodicity of the lattice, we look for solutions of (2) in terms of Floquet-Bloch waves (Floquet, 1883; Bloch, 1928; Bensoussan et al., 1978; Brillouin, 1978; Wilcox, 1978; Kittel, 1986). So, for a cubic array of unit cells Y,

\begin{array}{l} p (x + 1, y + 1, z + 1) = p (x, y, z) e^{i (k_{x} + k_{y} + k_{z})} & (3) \end{array}

where k_x, k_y and k_z are components of the Bloch vector k within the Brillouin zone Y* = [0, π]³ (Joannopoulos et al., 1995; Gazalet et al., 2013).

The implementation in the finite element package is fairly straightforward (Hladky-Hennion et al., 1991; Nicolet et al., 2004). We first multiply equation (2) by a smooth test function v and using integration by parts, we obtain the so-called weak form of the time-harmonic acoustic equation

\begin{matrix} - \int_{Y} ρ^{- 1} \nabla p \cdot \nabla v d x d y d z + \int_{\partial Y} ρ^{- 1} \frac{\partial p}{\partial n} v d s + ω^{2} \int_{Y} K^{- 1} p v d x d y d z = 0 & (4) \end{matrix}

where ∂f/∂n = ∇f·n with n the unit outward normal to the boundary ∂Y of Y, and ds the infinitesimal surface element on ∂Y.

We note that the weak formulation holds for heterogeneous fluids as ρ and K can be spatially varying. In particular, this model works for domains such as a homogeneous fluid filled with a periodic array of rigid obstacles. For the finite element implementation, (4) is discretised using test functions taking values on nodes of a tetrahedral mesh of the basic cell (first order tetrahedral elements), see e.g., (Nicolet et al., 2004) for further details. From (4), we note that setting rigid boundary conditions on an inclusion amounts to assuming Neumann (natural) homogeneous data, whereas transmission conditions at the interface between various fluid phases mean that the quantity ρ⁻¹∂p/∂n is preserved across the interface. We note that in the case of airborne pressure waves, the contrast in density between air and inclusions made of metal or even polymer is sufficiently large to assume Neumann data, but this simplification does not hold if we replace air by water.

Let us now consider a periodic array of defects Ω₁, …, Ω_N embedded in an elementary cell Y = ]0;1[³. Let p(x, y, z) satisfy the Helmholtz equation in $Y \ ⋃_{j = 1}^{N} \bar{Ω_{j}}$ . We also assume that p satisfies Neumann boundary conditions on the contours of defects, where n denotes the unit outward normal to the boundary ∂Ω_j of a defect Ω_j:

\begin{array}{l} \frac{\partial p}{\partial n} |_{\partial Ω_{j}} = 0, j = 1, \dots, N & (5) \end{array}

We would like to consider a particular case when the defects Ω₁, …, Ω_N are spherical shells with thin air channels connecting an air-filled interior cavity to the exterior surrounding air. These defects can be modeled as multistructures (Kozlov et al., 1999) in the following way,

\begin{array}{l} Ω_{(N)} = {a_{(N)} < \sqrt{x^{2} + y^{2} + z^{2}} < b_{(N)}} \ \bar{⋃_{j = 1}^{N} Π_{ε (N)}^{(j)}} & (6) \end{array}

where a_(N) and b_(N) are given constants and $Π_{ε (N)}^{(j)}$ is the thin channel.

4. Asymptotic Approximation: A Discrete Spring-Mass Model

In this section, we derive an asymptotic approximation of the field in thin channels $Π_{ε}^{(j)}$ , see Figure 3,

\begin{array}{l} Π_{ε}^{(j)} = {(x, y, z) \in R^{3} : 0 < x < l_{j}, \\ \sqrt{y {(t)}^{2} + z {(t)}^{2}} < ε h_{j} (t), (0 \leq t \leq 2 π)} & (7) \end{array}

where l_j is the length of the j^th bridge, εh_j(t) the radius of its varying cross-section D_ε (parametrized by a positive real t). Here, ε is a small positive non-dimensionnal parameter. To derive the asymptotic expansions, we introduce the scaled variables ξ = (y/ε, z/ε).

FIGURE 3

Figure 3. Geometry of the inclusions and the Helmholtz oscillator consisting of one spring connected to a mass at one end and fixed at the other end.

Without loss of generality, and for the sake of simplicity, we drop the superscript j. In Π_ε, the time-harmonic wave equation takes the rescaled form

\begin{array}{l} {\frac{1}{ρ} (\frac{1}{ε^{2}} Δ_{ξ} + \frac{\partial^{2}}{\partial x^{2}}) + \frac{ω^{2}}{K}} p = 0, & (8) \end{array}

with the Neumann boundary conditions

\begin{array}{l} \frac{\partial p}{\partial n} |_{\partial D_{ξ}} = 0 & (9) \end{array}

The field p is approximated in the form

\begin{array}{l} p ~ p^{(0)} (x, y, z) + ε^{2} p^{(1)} (x, y, z) & (10) \end{array}

To leading order, we obtain

\begin{array}{l} {\begin{array}{l} Δ_{ξ} p^{(0)} & = & 0 & on D_{ξ} \\ \nabla_{ξ} p^{(0)} & = & 0 & on \partial_{D_{ξ}} \end{array} & (11) \end{array}

Hence, p⁽⁰⁾ = p⁽⁰⁾(x) (it is ξ−independent). Assuming that p⁽⁰⁾ is given, we derive that the function p⁽¹⁾ satisfies the following model problem on the scaled cross-section of Π_ε

\begin{array}{l} {\begin{array}{l} Δ_{ξ} p^{(1)} & = & - \frac{1}{ρ} \frac{\partial^{2} p^{(0)}}{\partial x^{2}} + \frac{ω^{2}}{K} p^{(1)} & in D_{ξ} \\ \nabla_{ξ} p^{(1)} \cdot n & = & 0 & on \partial_{D_{ξ}} \end{array} & (12) \end{array}

The condition of solvability for the problem has the form:

\begin{array}{l} \frac{1}{ρ} \frac{d^{2} p^{(0)}}{d x^{2}} + \frac{ω^{2}}{K} p^{(0)} = 0, 0 < x < l_{j} & (13) \end{array}

Hence, we have shown that to the leading order we can approximate the field p in the thin channel Π_ε by the function p⁽⁰⁾ which satisfies the Helmholtz's equation in one-space dimension. We now assume that the field is periodic over the cell since it is localized. This shows that the average of the eigenfield over the macro-cell vanishes. Indeed, let χ₁ denote the value of the field in the large body Σ of the multi-structure Ω and let χ₂ (which we normalize to 1) denote the value of the field within the complementary area of the macro-cell Y \ Ω excluding the thin channels. Taking v = 1 in (4), we deduce that

\begin{array}{l} ω^{2} \int_{Y} ρ p d x d y d z = - \int_{\partial Y \cup \partial Ω} K \frac{\partial p}{\partial n} d S = 0 & (14) \end{array}

This shows that the average of the field p over Y vanishes, hence by neglecting the small volume of the thin channels, we obtain

\begin{array}{l} χ_{1} {meas}_{Σ} + χ_{2} {meas}_{Y \ Ω} = O (ε) & (15) \end{array}

where meas_Σ and meas_{Y \ Ω} denote, respectively the volumes of Σ and Y \ Ω.

We now consider two cases. The first one is the study of an array of simple spherical shells with either one or six thin channels, and the other one is the study of an array of double spherical shells with one thin channel in each shell. Since we have q thin channels, we have q separate eigensolutions V_j, (j = 1, …, q), corresponding to the vibrations of thin domains $Π_{ε}^{(j)}$

\begin{array}{l} ρ^{- 1} V_{j}^{″} (x) + K^{- 1} ω^{2} V_{j} (x) = 0, 0 < x < l_{j}, & (16) \end{array}

\begin{array}{l} V_{j} (0) = χ_{2} = - χ_{1} \frac{meas (Ξ)}{meas (Y \ Ω)} . & (17) \end{array}

We note that V_j(0) is equal to a non-zero constant and also that V_j(l_j) = χ₁. Next, we need to take into account that some boundary layers occur at the end regions of thin ligaments. These boundary layers are characterized by exponential decay when the boundary conditions for the functions V₀ (the leading term in the asymptotic expansion of the thin bridge solutions) are chosen in a specific way. In our case, integrating (2) over $Σ \cup Π_{ε}^{(1)} \cup \dots \cup Π_{ε}^{(q)}$ and applying the divergence (or Gauss) theorem, we obtain to order O(ε)

\begin{array}{l} K^{- 1} I_{j} V_{j}^{'} (l_{j}) = M_{j} ω^{2} V_{j} (l_{j}), & (18) \end{array}

where

\begin{array}{l} I_{j} = \int_{0}^{2 π} ε h_{j} (t) d t . & (19) \end{array}

All the channels are connected to Ξ, hence, V₁(l₁) = … = V_q(l_q) = V. We note that the boundary layer condition (18) can be interpreted as Newton's second law.

The solution of the problem (16)−(18) (that one can interpret in terms of a spring-mass model like in Figure 3) has the form

\begin{array}{l} V_{j} (x) = - \frac{χ_{2} [cos ((ω / c) l_{j}) - 1]}{sin ((ω / c) l_{j})} sin (\frac{ω}{c} x) + χ_{2} cos (\frac{ω}{c} x) & (20) \end{array}

where $c = \sqrt{K / ρ}$ and the frequency ω is given as the solution of the following equation:

\begin{array}{l} \sum_{j = 1}^{n} (I_{j} cot (\frac{ω l_{j}}{c})) = \frac{m_{j} c}{K} ω & (21) \end{array}

where we invoked (18). Looking at a first low frequency, we deduce an explicit asymptotic approximation

\begin{array}{l} ω ~ \sqrt{\sum_{j = 1}^{n} (\frac{I_{j}}{l_{j}})} \sqrt{\frac{K}{M}} (1 + \frac{meas (Ξ)}{meas (Y \ Ω)}) & (22) \end{array}

This estimate actually holds for the frequency ω₂ of the upper edge of the phononic band gap. We note that if we take V(0) = 0 instead of V(0) = meas(Ξ)/meas(Y \ Ω), we estimate the frequency of the lower-edge of the phononic band gap. We also notice that the boundary layer condition (18) is only valid at order 0(ε), so it needs to be refined if one wants to improve the frequency estimate in (22).

4.1. Eigenfrequency Estimate in the Case of a Single Spherical Shell With one or six Thin Channels

We report in Figures 4, 5 finite element computations for a periodic cell of sidelength d = 1 m with a simple spherical shell with one thin channel. We then proceed with the same shell with six thin channels. The geometry of the elementary cell and associated spring-mass model are given in Figure 6 and finite element computations are shown in Figures 7, 8. The interior and exterior radii of the shell are respectively 0.3 and 0.4 m, the thin channels have the same length 0.1m and radii 0.01 m. Therefore, the frequency estimates are (in rad.s⁻¹):

\begin{array}{l} ω_{1} ~ 52.1107, ω_{2} ~ 54.0796 & (23) \end{array}

for one thin channel, which are in good agreement with the finite element values

\begin{array}{l} ω_{1}^{*} = 51.4250, ω_{2}^{*} = 55.2638 & (24) \end{array}

for one thin channel, and

\begin{array}{l} ω_{1} ~ 127.6447, ω_{2} ~ 132.4676 & (25) \end{array}

for six thin channels, which are in good agreement with the finite element values

\begin{array}{l} ω_{1}^{*} = 125.6586, ω_{2}^{*} = 134.8155 & (26) \end{array}

for six thin channels.

FIGURE 4

Figure 4. Dispersion diagram for a periodic array (pitch d = 1m) of spherical rigid shells (inner radius 0.3 m and outer radius 0.4 m) with one thin channel (length 0.1 m and radius 0.01 m) representing the frequency ω (rad.s⁻¹) of pressure waves in air vs. the wavenumber k = |k| (m⁻¹), projection of the Bloch vector k along the edges of the irreducible Brillouin zone ΓXMU shown in Figure 1. We note the appearance of a frequency stop band for ω ∈ [51.4250, 55.2638] rad.s⁻¹.

FIGURE 5

Figure 5. The eigenfunction corresponding to the eigenfrequency $ω_{1}^{*} = 51.4250$ rad.s⁻¹ for one thin channel. Blue color corresponds to nearly vanishing amplitude of the eigenmode u, while red color corresponds to it maximum value. The pressure field p is constant inside the inner cavity and outside the shell, but it varies rapidly inside the thin channel: it is a localized eigenmode responsible for the stop band in Figure 4, which is well approximated by a spring mass model.

FIGURE 6

Figure 6. Geometry of the inclusions and the Helmholtz oscillator consisting of six springs connected to a mass at one end and fixed at the other end.

FIGURE 7

Figure 7. Dispersion diagram for a periodic array (pitch d = 1 m) of spherical rigid shells (inner radius R = 0.3 m and outer radius R = 0.4 m) with six thin channels (length 0.1 m and radius 0.01m) representing the frequency ω (rad.s⁻¹) of pressure waves in air vs. the wavenumber |k| (m⁻¹), projection of the Bloch vector k along the edges of the irreducible Brillouin zone ΓXMU. We note the appearance of a frequency stop band for ω ∈ [125.6586, 134.8155] rad.s⁻¹ which is wider and at higher frequencies than the stop band in Figure 4: the more identical thin channels, the higher the resonant frequency of the localized mode.

FIGURE 8

Figure 8. The eigenfunction corresponding to the eigenfrequency $ω_{1}^{*} = 125.6586$ rad.s⁻¹ for six thin channels responsible for the stop band in Figure 7. This frequency is well approximated by the spring mass discrete model which provides us with the frequency estimate ω₁ = 127.6447 rad.s⁻¹.

This demonstrates that the discrete model provides accurate estimates for the lower and upper edges of the ultra-low frequency stop band. This is therefore a useful tool which can be used in the design of acoustic metamaterials. Importantly, we note that the spring-mass counterpart of the shell with 6 holes, is like in Figure 6, which corresponds to an LC electrical circuit with one capacitance C and six inductances L in series corresponding to the mass and springs, respectively, see e.g., Guenneau et al. (2007a) for an analysis of two-dimensional acoustic metamaterials.

4.2. Dependence of low Frequency Stop Band on Radius of Thin Channel

Before we move to more complex geometries, we would like to investigate the case of Helmholtz resonators with thicker channels, that would prevent viscous effects for pressure waves entering the thin channels in manufactured prototypes. We note that there is an obvious limitation of the asymptotic model. Indeed, the accuracy of the eigenfrequency estimate relies heavily on the smallness of the parameter ϵ, which essentially requires elongated channels (in other terms ligaments or necks). We have checked the lack of robustness of the asymptotic estimates with respect to the elongation of channels, when we changed the radius of the thin channel for a Helmholtz resonator with a single channel like in Figure 3. We therefore focus on finite element results which are shown in Table 1, and one can see that although the low frequency stop band is preserved for increasing values of the radius, its position is shifted up and it is enlarged. These results are good news for forthcoming experiments.

TABLE 1

Table 1. Values of the extrema of the gap with respect to the radius of the thin channel, for the Helmholtz resonator like in Figure 3.

FIGURE 9

Figure 9. Modification of the gap width and position with respect to the radius of the thin channel of the Helmholtz resonator with a single shell like in Figure 3. The eigenfrequency estimates ω₁~52.1107 and ω₂~54.0796 compare well with finite element results $ω_{1}^{*} ~ 51.4250$ and $ω_{2}^{*} ~ 55.2638$ for r = 0.01m, see Table 1. Note the non linear scale for r.

4.3. Eigenfrequency Estimate in the Case of a Double Spherical Shell With One Thin Channel

Let us now consider the case of and LC circuit with more than one capacitance. The simplest model of interest is that of an LC circuit with two capacitances and two inductors. This corresponds to a mechanical model like in Figure 10.

FIGURE 10

Figure 10. Geometry of the inclusions and the Helmholtz oscillator consisting of two masses connected by a spring, with one of them connected to a fixed domain.

In the numerical example we now have $ε^{2} h_{2}^{2} = 3.14 \times 1 0^{- 4}$ m², $ε^{2} h_{1}^{2} = 7.85 \times 1 0^{- 5}$ m², l₂ = 0.1m, l₁ = 0.05m, and the masses (in kilogram)

\begin{array}{l} m_{1} = ρ V_{1} = \frac{4000}{3} π r_{1}^{3} \\ m_{2} = ρ V_{2} = \frac{4000}{3} π (r_{1}^{3} + (b_{2}^{3} - a_{2}^{3})) + 1 0^{3} ε^{2} h_{1}^{2} l_{1} & (27) \end{array}

where ρ is the density of air (~1.225kg.m⁻³), V₁ and V₂ the volumes air occupies in Ξ₁ and Ξ₂, r₁ is the interior radius for the domain Ξ₁ and a₂, b₂ are respectively the interior and exterior radii for the domain Ξ₂. In our case, r₁ = 0.1m, a₂ = 0.15m and b₂ = 0.2m. The formula (22) gives the following values for the first eigenfrequencies (in rad.s⁻¹) of the multistructures Π_ε(1)⋃Ξ₍₁₎ and Π_ε(2)⋃Ξ₍₂₎:

\begin{array}{l} ω_{1} ~ 124.2641, ω_{2} ~ 208.4109 & (28) \end{array}

The corresponding angular frequencies (in rad.s⁻¹) associated with the standing waves in the periodic structure were obtained numerically, and from the band diagram in Figure 11 they are

\begin{array}{l} ω_{1}^{*} = 55.2128, ω_{2}^{*} = 118.9750 & (29) \end{array}

FIGURE 11

Figure 11. Dispersion diagram for a periodic array (pitch d = 1m) of double spherical rigid shells (radius of spheres from inner to outer are 0.1, 0.15, 0.3, and 0.4 m) with one thin channel in each shell (respectively of lengths 0.05 and 0.1 m and radii 0.005 and 0.01 m) representing the frequency ω (rad.s⁻¹) of pressure waves in air vs. the wavenumber k = |k| (m⁻¹), projection of the Bloch vector k along the edges of the irreducible ΓXMU. We note the appearance of two frequency stop bands for ω ∈ [55.2128, 58.5863] rad.s⁻¹ and ω ∈ [118.9750, 119.0959] rad.s⁻¹.

Formula (22) leads to frequency estimates in (28) that do not capture the eigenfrequency $ω_{1}^{*}$ , that corresponds to the localized eigenmode shown in the left part of Figure 12, but we observe a good agreement between frequency estimate ω₁ and $ω_{2}^{*}$ , associated with the eigenmode shown in right part of Figure 12. The estimate for the eigenfrequency $ω_{1}^{*}$ can be found if the domain Π_ε(2)⋃Ξ₍₂₎ is replaced by the domain Π_ε(2)⋃Ω₍₂₎⋃Π_ε(1)⋃Ω₍₁₎. In this case, the eigenfrequency $ω_{1}^{*}$ is approximated by the first positive eigenvalue of the problem

\begin{array}{l} ρ^{- 1} V_{1}^{″} (x) + K^{- 1} ω^{2} V_{1} (x) = 0, 0 < x < l_{1} & (30) \end{array}

\begin{array}{l} V_{1} (0) = 0, & (31) \end{array}

\begin{array}{l} K^{- 1} I_{1} V_{1}^{'} (l_{1}) - K^{- 1} I_{2} V_{2}^{'} (0) = m_{1} ω^{2} V_{1} (l_{1}), & (32) \end{array}

\begin{array}{l} ρ^{- 1} V_{2}^{″} (x) + K^{- 1} ω^{2} V_{2} (x) = 0, 0 < x < l_{2} & (33) \end{array}

\begin{array}{l} λ^{- 1} I_{2} V_{2}^{'} (l_{2}) = m_{2} ω^{2} V_{2} (l_{2}), & (34) \end{array}

\begin{array}{l} V_{2} (0) = V_{1} (l_{1}), & (35) \end{array}

where V₁(x), V₂(x) are the eigenfunctions defined on (0, l₁) and (0, l₂), respectively, and the masses m₁, m₂ are defined by (in kilogram)

\begin{array}{l} m_{1} = \frac{4000}{3} π r_{1}^{3}, m_{2} = \frac{4000}{3} π (b_{2}^{3} - a_{2}^{3}) & (36) \end{array}

FIGURE 12

Figure 12. Eigenfunctions for embedded resonators corresponding to the eigenfrequency $ω_{1}^{*} = 55.2128$ rad.s⁻¹ (left panel) and corresponding to the eigenfrequency $ω_{2}^{*} = 118.975$ rad.s⁻¹ (right panel). In left panel both canals vibrates, while in right panel only the inner canal does. The second eigenfrequency $ω_{2}^{*} = 118.975$ rad.s⁻¹ is reasonably well approximated by the first eigenfrequency of the spring mass model which gives and ω₁ = 124.2641 rad.s⁻¹ and ω₂ = 208.4109 rad.s⁻¹.

Taking into account that ω₀ = O(ε), we deduce that it can be approximated as the first positive solution of the following algebraic equation:

\begin{matrix} m_{1} m_{2} l_{1} l_{2} ω^{4} - K ω^{2} (l_{2} I_{1} m_{1} + I_{1} l_{2} m_{2} + I_{2} m_{1} l_{1}) + K^{2} I_{1} I_{2} = 0 & (37) \end{matrix}

so that ω₀~82.76057 rad.s⁻¹, which provides a reasonably accurate approximation of ω₁* = 55.2128 rad.s⁻¹. However, a further refinement of this asymptotic estimate would require adding higher order corrections and thus solving boundary layer type problems as discussed in Kozlov et al. (1999), and this falls beyond the scope of the present paper.

5. Isofrequency Surfaces, Dynamic Effective Anisotropy and Localized Modes

Thus far, we have only analyzed dispersion curves, but it has been known for over 80 years that one needs to be extra careful regarding analysis of say electronic band structures (Bouckaert et al., 1936), as there are examples showing that stop band edges might be reached strictly within the Brillouin zone, and so it is not enough to describe its edges to characterize the band spectrum (Harrison et al., 2007). We would like to investigate dynamic effective anisotropy effects near the edge of the low frequency stop band, and localization effects induced by local resonators such as in Movchan et al. (2006) ,Bigoni et al. (2013) ,Craster et al. (2013), and Llewellyn Smith and Davis (2010) for a frequency on the acoustic band.

Here, we provide representation for isofrequencies near the lower edge of the first stop band for the case of a Helmholtz resonator with one thin channel as in Figure 4. However, we observed very similar features in all other cases. We first note in Figure 13 that the distortion of the isofrequency surfaces increases when we move toward the lower edge of the stop band. A small change in the frequency leads to a dramatic change in the isofrequency surface. Besides, the observation of isofrequency (or slowness) contours in Figure 13 confirms that anisotropy increases near the edge of the stop band as contours become more and more elongated. Researchers in photonics used vanishing group velocity near stop band edges to achieve self-guiding of light in 2D photonic crystals in the past 20 years (Kosaka et al., 1999; Witzens et al., 2002; Chigrin et al., 2003). Such an effect is known in the applied mathematics and waves literature as a dynamic effective anisotropy (Slepyan et al., 1987; Ayzenberg-Stepanenko and Slepyan, 2008; Craster et al., 2009, 2010; Colquitt et al., 2012; Schnitzer, 2017) : for an isotropic medium, the isofrequency surface is spherical, and the corresponding isofrequency contours are circular. On the edge of the stop band, the extreme elongation of isofrequency contours reveals an extreme dynamic anisotropy. We note that some dynamic effective anisotropy of 3D dynamic lattices displays similar isofrquency surfaces to ours near stop band edges (Vanel et al., 2016). This could be used to achieve as aforementioned highly directive phenomena, in a way similar to what was demonstrated experimentally for transverse electromagnetic waves (Ceresoli et al., 2015). Another striking effect is that of wave localization, when one creates one or more defects in a periodic structure. Usually, this is achieved in high frequency stop bands, when the wave wavelength is on the same order as the periodicity. Here, the wave localization is achieved either at resonances of resonators, see upper panel in Figure 14, or as shown in the lower panel of the same figure, at a frequency $ω_{1}^{*} = 559.2042$ rad.s⁻¹ located in a part of the acoustic band sandwiched between the low frequency stop band and another region of vanishing group velocity near the X symmetry point in Figure 4, so this has the additional feature of being below the Bragg frequency regime. We note that from Figure 14, this localized mode can be approximated as a monopole given by the first derivative of the spherical Bessel function of the first kind j_n(ωr/c), where r is the radial position. The lowest root corresponding to the lowest frequency of the cavity is the frequency given by the root of the transcendental equation (simply derived by separation of variables in Helmholtz's equation (2) written in spherical coordinates and assuming that ∂p/∂n = 0 on the boundary of the defect) :

\begin{array}{l} j_{1}^{'} (ω a / c) = 0 . & (38) \end{array}

Bearing in mind that the first root of $j_{1}^{'}$ is 2.0816 and taking a = 1.25m as the radius of a dashed circle within which the defect mode seems to be mostly confined in Figure 14, (38) leads to the frequency estimate ω₁ = 566.19 rad.s⁻¹, which is in good agreement with the finite element computation $ω_{1}^{*} = 559.2042$ rad.s⁻¹. This method of the effective sphere is the analogous concept to effective disk in doubly periodic structures with defects (Movchan et al., 2007).

FIGURE 13

Figure 13. Representative isofrequency surfaces k_z(k_x, k_y) in the first Brillouin zone (upper panel) and associated isofrequency contours in the plane (k_x, k_y) for k_z ∈ (0, π) (lower panel) at frequency ω = 51.42 rad.s⁻¹ (left panel), ω = 51.4215 rad.s⁻¹ (middle panel) and ω = 51.423 rad.s⁻¹ (right panel); The surfaces and contours flatten when the frequency moves toward the lower edge of the stop band in Figure 4. and thus dynamic effective anisotropy increases.

FIGURE 14

Figure 14. Representative eigenfunctions for a macrocell of 26 resonators as in Figures 4, 5 with a defect (air instead of resonator) in the middle. Upper panel: Eigenfunctions corresponding to an eigenfrequency $ω_{0}^{*} = 52.3545$ rad.s⁻¹ in the ultra-low frequency stop band of Figure 4. Lower panel: Eigenfunction corresponding to an eigenfrequency $ω_{1}^{*} = 559.2042$ rad.s⁻¹ above the low frequency stop band of Figure 4, which is approximated by the first root of (38) i.e., the fundamental resonance of an effective spherical cavity (marked by dashed circle) of radius a = 1.25 m.

6. Conclusion

In this paper, we have seen that it is possible to sculpt the Bloch spectrum of three-dimensional phononic crystals almost ad libitum simply by digging some holes and adding cavities in rigid spheres periodically arranged along a cubic lattice. One of the main achievements of our numerical study is the appearance of ultra-low frequency stop bands at frequencies predicted quantitalively by an asymptotic model, that can be viewed as a 3D counterpart of (Movchan and Guenneau, 2004). We also conducted some elementary shape optimization (by varying the size, diameter and number of channels in a rigid sphere of constant radius) in order to enhance the control of the location and the number of low frequency stop bands, thanks to our asymptotic estimates. Our asymptotic results could be further used for a homogenization study in the spirit of the 2D analysis conducted in Mercier et al. (2017) to address the effective parameters of our acoustic system. However, the acoustic metamaterial we have analyzed would allow only frequency dependent effective bulk modulus, and in order to achieve a frequency dependent effective density and thus double negative acoustic parameters near resonances, one needs to introduce a second type of resonator in the periodic cell. Therefore, achieving a design of an acoustic superlens through negative refraction requires a more complex design. However, a cubic array of rigid spheres does not support any complete stop band, even in the densely packed configuration, so the acoustic metamaterial which we studied has markedly different filtering properties compared with phononic crystals with rigid spheres: multistructures make possible tunable ultra-low frequency stop bands (associated with very flat dispersion curves i.e., localized eigenmodes). Finally, we illustrate in Figure 14 a possible application of the ultra-low frequency stop band in order to localize a mode of a wavelength larger than the pitch of the array of resonators. One can also envisage sculpting a line defect thanks to the removal of a few resonators in a larger macrocell, that would make a low frequency waveguide. We further note that considering shells of same diameter with holes of varying diameters, or spheres with varying diameters and identical or varying holes would make possible graded phononic crystals with a larger low frequency stop band, thanks to the overlap of numerous low frequency stop bands. This could lead to a three dimensional counterpart to the acoustic rainbow put forward in Zhu et al. (2013). Indeed, many embedded spheres would make this possible, since the band diagram 11 suggests N embedded spheres would lead to N low frequency stop bands that might hybridize. Similarly, one could envisage to reflect, detour, or focus, pressure waves using the low frequency stop band within which effective parameters are expected to take negative values as it is now well-established for the two-dimensional counterpart of such kind of acoustic metamaterials. We hope our work will foster experimental efforts toward 3D acoustic metamaterials for airborne and underwater sound filtering effects. At first sight, the implementation of the former seems more straightforward, since the latter would require further theoretical and numerical investigation of conversion of pressure waves propagating in the fluid into pressure and shear waves propagating in the shells. For instance, in the derivation of the asymptotic estimates, the boundary layer analysis will be more involved. Nevertheless, it might be easier to perform underwater experimental characterization of the sonic metamaterial.

Author Contributions

GD performed all numerics, AM derived the asymptotic formulae, SG wrote the manuscript and SE participated in the physical interpretation of results. All authors reviewed the manuscript and contributed to the study.

Funding

This research has been funded in part by ERC (SG received funding from ERC as a PI of an ERC starting grant 279673 during the period 2011–2016, with GD a postdoc on that grant in 2014-2015) and by EPSRC (AM received funding from EPSRC as a co-PI of research grant EP/L024926/1 during the period 2015-2018 and SG is a named collaborator who is currently hosted in the group of Prof. R. Craster at Imperial College who is the PI on that grant).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

AM acknowledges funding of the EPSRC as a co-investigator on program grant EP/L024926/1 on Mathematical fundamentals of Metamaterials for multiscale Physics and Mechanics. SG is thankful for a visiting position in the department of mathematics at Imperial College London and support from EPSRC as a named collaborator on grant EP/L024926/1. Insightful discussions with Prof. R.V. Craster are acknowledged.

References

Aberg, M., and Gudmundson, P. (1997). The usage of standard finite element codes for computation of dispersion relations in materials with periodic microstructure. J. Acoust. Soc. Am. 102, 2007–2013.

Google Scholar

Alster, M. (1972). Improved calculation of resonant frequencies of helmholtz resonators. J. Sound Vib. 24, 63–85.

Google Scholar

Ammari, H., Fitzpatrick, B., Lee, H., Yu, S., and Zhang, H. (2017b). Double-negative acoustic metamaterials. arxiv:1709.08177v2.

Google Scholar

Ammari, H., and Zhang, H. (2017a). Effective medium theory for acoustic waves in bubbly fluids near Minnaert resonant frequency. SIAM J. Math. Anal. 49, 3252–3276. doi: 10.1137/16M1078574

ORIGINAL RESEARCH article

Analysis of Low Frequency Acoustic Stop Bands in Cubic Arrays of Thick Spherical Shells With Holes

1. Introduction: Acoustic Metamaterials

2. Motivation: Spectral properties of a periodic array of rigid spheres

3. Setup of the spectral problem: The continuum model

4. Asymptotic Approximation: A Discrete Spring-Mass Model

4.1. Eigenfrequency Estimate in the Case of a Single Spherical Shell With one or six Thin Channels

4.2. Dependence of low Frequency Stop Band on Radius of Thin Channel

4.3. Eigenfrequency Estimate in the Case of a Double Spherical Shell With One Thin Channel

5. Isofrequency Surfaces, Dynamic Effective Anisotropy and Localized Modes

6. Conclusion

Author Contributions

Funding

Conflict of Interest Statement

Acknowledgments

References

This article is part of the Research Topic

People also looked at