Uncertain acoustic meta-atoms
Felix Kronowetter(1) , Martin Eser(1) , Kheirollah Sepahvand(1) , Steffen Marburg(1)(1) Technical University of Munich, Chair of Vibroacoustics of Vehicles and Machines, Germany, felix.kronowetter@tum.de

Abstract
Acoustic metamaterials (AMMs) consist of periodic arrangements of single meta-atoms (e.g. Brillouin zones).
The AMMs can manipulate the acoustic wave propagation in ways that are not found in nature or conventional
materials. Furthermore, AMMs can have unnatural material properties such as a negative effective mass or bandgaps. One type of meta-atom is based on the principle of a Helmholtz resonator that is embedded in a fluid
matrix. A periodic arrangement of such meta-atoms in the two dimensional space combines the effects of a resonator and those of phononic crystals. The effectiveness of that kind of AMM depends on the eigenfrequenciesof the resonators and the relative position of one meta-atom to one another. Since the production of AMMs
is linked to manufacturing tolerances the perfect periodicity is not fulfilled and can affect the properties of theAMM. This work deals with the uncertainties of the meta-atoms concerning the geometry of the embedded
resonator. The uncertain geometry parameters are approximated by spectral expansions combined with the nonintrusive collocation method. Further, the transfer function and the insertion loss with respect to the uncertain
parameters are analyzed. Finally, the results of the spectral approach are compared to those of the Monte Carlo
method.
Keywords: Acoustic metamaterials, Uncertainty quantification, Infinite element method, Spectral expansions

1

INTRODUCTION

The demand for acoustic noise barriers has increased over the last years due to multiple environmental regulations like the reduction of the drive by noise of cars or the need of suppressing the propagation of disturbing
noise frequencies in airplanes. Sonic crystals and their corresponding acoustic Bragg type band gaps have encountered a lot of attention when it comes to the design of noise barriers [1],[2],[3]. Since the band gaps ofsonic crystals are strongly dependent on the lattice constant of the periodic arrangement [4] they can be extended by resonant elements [4],[5]. Hence, the structure of resonant elements does not only feature the Bragg
type band gaps, but also additional ones, that can be assigned to the resonant frequencies of the resonators.
Several solution techniques exist for the evaluation of the resonant frequencies. One is based on the assumption, that the structure consists of infinite arrays of cylinders or resonators [5]. Accordingly, the whole problem
is transferred onto a reciprocal lattice and a unit cell is defined. Based on the unit cell the Floquet-Blochtheorem can be applied, that results in periodic boundary conditions depending on the wave vector in the kspace, that is defined in the first Brillouin zone. The evaluation of eigenfrequencies using the Floquet-Blochtheorem yields the so called dispersion curves that can visualize the band gaps in the k-space. This techniqueneglects the finite dimensions of the array in a real application and also the effects occurring at the outer edgesof the structure. Another solution technique implicates that the examined array of cylinders or resonators ismuch smaller than their surroundings. Therefore, the problem is denoted an unbounded or exterior Helmholtz
problem. The encircling fluid domain is assumed to be of infinite extension and the non-reflective Sommerfeld
radiation condition [6] is applied to the outer boundary. The unboundedness of the problem leads to difficulties
in numerical solution techniques. To overcome these difficulties, the infinite element method [7] is applied in
combination with the finite element method. This combined method provides the mass, damping and stiffness
matrices to compute the quadratic eigenvalue problem. The identification of the resonant modes, also cavity or
trapped modes, is done by the quality factor of the complex eigenvalues. Another aspects considered in thispaper are manufacturing tolerances, that result in a derivation of the real geometry compared to the computer

1963model. The manufacturing tolerances are represented by independent random variables. That is why the infinite element method is extended by spectral expansions. The random input data lead to somehow distributed
eigenvalues / resonant frequencies. Since the distribution of the resonant frequencies is dependent on the random input data, the probability density function defines the possible frequency range. In this paper just two
dimensional problems and single c-shaped resonators, also called meta-atoms, are examined. The procedure issimilarly applicable for periodic arrangements of the meta-atoms, furthermore for three dimensional problems.

2

MATHEMATICAL AND GEOMETRICAL MODELS

2.1 Infinite element method
The effects described by linear acoustics in a fluid domain Ω f are covered by the elliptic Helmholtz-equation(1) in the frequency domain. In case of a time-harmonic excitation, the steady solution is formulated as p̃(xx,t) =ℜ{p(xx)e−iωt }, with p̃(xx,t) being the sound pressure. The time-harmonic sound pressure p is denoted p = p for
convenience.
∆p(xx) + k2 p(xx) = 0 , x ∈ Ω f ⊂ Rd
(1)
The wavenumber is denoted k = ωc and the circuar frequency ω = 2π f .
For interior Helmholtz problems, the finite element discretization of equation (1) applying the Bubnov-Galerkinmethod and the Robin boundary condition leads to the following system of equations in matrix form [8]K − ikC
C − k2M )pp = f(K(2)

with frequency independent matrices. For exterior Helmholtz problems the fluid domain is theoretically of
infinite extension. To be able to solve the problem numerically the fluid domain is separated into a boundeddomain and a complementary domain Ω fc . This is exemplarily shown in figure 1.
Ω fc
Ωf
Γc

Γ∞

Figure 1. Bounded fluid domain Ω f closed by Γc and unbounded fluid domain Ω fc with the boundary Γ∞ atinfinity.
In the complementary domain Ω fc the infinite element method (IFEM) is applied. This method also solves theHelmholtz equation and furthermore fulfills the Sommerfeld radiation condition [6]. For the infinite elements,
the conjugated Astley-Leis formulation is chosen. This formulation is based on the Petrov-Galerkin method,
does not require a truncation of the domain due to the mapping of the infinite elements and leads to frequency
independent matrices [7]. The choice of the weight and test functions being complex conjugates results in amore accurate resolution of the complementary domain compared to the unconjugated formulation [9],[10],[11].
The coupled FE-IFE system of equation reads as:
K f fc − ikC
C f fc − k2M f fc )pp = f .
(K

1964

(3)


Where K f fc =Kf
K fc fK f fc
K fc




,

C f fc =Cf
C fc fC f fc
C fc




,

M f fc =Mf
M fc f

M f fc.
M fc

The matrices with the indices f fc and fc f contain the entries, that result from the coupling nodes at theboundary Γc between the bounded fluid domain and the complementary domain.
2.2 Uncertainty quantification
Most simulations use deterministic input data that correspond to a certain reference or mean value. In this paper
we assume the input data to be randomly distributed due to manufacturing tolerances. Hence, the input data isparameterized ξ = {ξ1 , . . . , ξN } with ξi (θ ) being vector quantities in the size of elements in the sample space
i = 1, . . . , Ns .
One solution method to problems with random input data is the Monte-Carlo method. This method leads to an
exact solution but the problem has to be solved for any realization of the input data. Thus, the method requires
extensive computational power. Little realizations of the input data result in a more inaccurate approximationof the mean value and of the variance of the solution. The mathematical convergence of the mean value of thesolution applying the Monte-Carlo method reads as ∼ √1M , where M is the number of realizations [12].
Spectral methods on the other side are more efficient and advanced however more complex in the application.
Spectral expansions are based on Fourier-like expansions of a random process in L2 , that are convergent withrespect to the norm associated with the corresponding inner product. The polynomial chaos expansion first
introduced by Wiener [13] is restricted to Hermite polynomials spanning the orthogonal base. The choice ofpolynomials for the expansion is dependent on the type of distribution of the random input data [14]. A moresophisticated model is the generalized polynomial chaos (gPC) [15], that is not limited to Hermite polynomials,
respectively to Gaussian distributed input data:∞X(ξξ ) = ∑ xi Φi (ξξ ).

(4)

i=0

The gPC is accurate related to the mean value of the data. For numerical calculations the infinite series ofequation (4) has to be truncated. The truncated, finite series (6) is dependent on the number of random variables
N and on the highest order of the polynomials p.
P+1 =(N + p)!
N!p!

(5)

P

X(ξξ ) = ∑ xi Φi (ξξ ) + ε(N, p)

(6)

i=0

The system of equations with spectral expansions of the random input data can be solved using either intrusive
or non-intrusive methods. In this paper solely non-intrusive methods are examined, where the system of equations can be treated as black-box. One type of non-intrusive methods is the so called collocation method, that
requires the system of equations to be solved at merely certain particularly chosen collocation points. Those
collocation points are defined by the roots of the polynomial of the order p + 1, are sorted referring to the region of the highest probability and the necessary amount of points determined [16]. The procedure is visualizedin equation (7). The deterministic spectral modes are denoted xi and the particular solutions of the system ofcpi ) for the corresponding collocation points cp i .
equations X ∗ (cp   ∗ 0 cp0 ) Φ1 (cpcp0 ) · · · ΦP (cpcp0 ) cp ) Φ0 (cpX (cp x0 





1)
1)
1)  
∗ (cp Φ0 (cp




cp
cp
cp
cp1 ) 
Φ
(cp
·
·
·
Φ
(cp
X
x
P
1

 1
=
(7)


..
..
.
.
.
.
..
..
..
.. 




.
.








∗
∗
∗
   ∗ P∗ xP
cpP ) Φ1 (cpcpP ) · · · ΦP (cp
cpP )
cp )
Φ0 (cpX (cp1965

The extension of equation (3) to the dependence of the matrices on random variables leads to:noC(·,ξξ ) − k2M (·,ξξ ) p (·,ξξ ) = f (·,ξξ ).
K (·,ξξ ) − ikC(8)

The application of spectral expansion to the system of equations (8) results in the spectral-stochastic infinite
element method (SSIFEM):
n T
o
T
T
T
Φ(ξξ ) − ikC
C (·)Φ
Φ(ξξ ) − k2M (·)ΦΦ(ξξ ) p T (·)ΦΦ(ξξ ) = f (·)ΦΦ(ξξ ).
K (·)Φ
(9)
2.3 Modal analysis and quality factor
The computation of normal modes can be done applying the finite element method and infinite element method[17]. Regarding equation (3) this leads to a quadratic eigenvalue problem (10).
M + λC̃C + K̃
K
Q (λ ) = λ 2M̃

(10)

K are the adapted matrices from equation (3). Both left and right eigenvectors are computed
C , K̃
M ,C̃
The matrices M̃
and the full problem is solved [18]. To be able to solve a quadratic eigenvalue problem with random input data
and applying the generalized polynomial chaos results in a spectral-stochastic formulation:
TTT

Φ(ξξ ) = λ 2M̃
M (·)Φ
Φ(ξξ ) + λC̃
C (·)Φ
Φ(ξξ ) + K̃K (·)Φ
Φ(ξξ ).
Q T (λ , ·)Φ

(11)

The eigenvalues to that kind of problem are complex numbers in the complex plane. Hence, the quality factor}of the eigenvalues can be computed by using q f ac ∝ | ℜ{λℑ{λ } | [19],[20].
2.4 C-shape meta-atom
The two dimensional c-shaped resonator studied by [4], [5] and [17] is used for the modal analysis. The exact
geometry is shown in figure 2.

r(ξ j )rA

A(ξi )

Figure 2. C-shaped resonator with with deterministic parameters (left) and random parameters (right).
For the two dimensional c-shaped resonator the outer radius r and the sloth width A are assumed to be randomly
distributed due to manufacturing tolerances. Both random variables are considered independent and standard
normal distributed. Therefore, Hermite polynomials are used in the expansions to model the input and output
data. The outer radius is set to r = 0.02m, the wall thickness to dr = 0.005m and the slot width to A = 0.005m.

1966

3

RESULTSThe eigenfrequencies of the uncertain c-shaped meta-atom are calculated using MATLAB. To safe computational power, the domain is formulated as half-space problem. The x-axis is chosen as symmetry axis and theadmittance boundary set to zero Y = 0 to create a fully reflecting, sound-hard wall. The boundary between thec-shaped resonator and the fluid domain is also assumed to be sound hard. This is shown in figure 3. For
the modal analysis quadratic elements are chosen for the finite element region and also quadratic interpolation
functions for the transversal direction of the infinite elements, whereas the radial order of the infinite elements
is 6 and Lagrange polynomials are used.

Ω fc

Ωf

Y =0
Figure 3. Half-space c-shaped meta-atom.
The sound pressure in the finite element region Ω f for the first 4 resonant frequencies, frequencies with anascending real part and a descending quality factor, of the meta-atom are shown in figure 4. As it can be seen,
the sound pressure is trapped in the meta-atom and shows little radiation to the far field. Hence, the quality
factor is a good criteria to find the resonant frequencies in the complex plane. The random input data is set
to r(ξ1 ) = rµ + σr ξ1 for the random radius and to A(ξ2 ) = Aµ + σA ξ2 , with rµ = r, σr = 0.05 · r, Aµ = A andσA = 0.1·A. The number of random variables is n = 2 and the order of the polynomial of the expansion is p = 3.
The random input data are chosen randomly, without any experimental background and are over dimensioned
to visualize their influence on the output data. Nevertheless, any configuration of the input data is possible.
The probability density functions (pdfs) for the random eigenfrequencies of the c-shaped resonator are shown infigure 5. The yellow line marks the mean value of the data, that is in accordance with the deterministic resultsof 4, meaning that the results are mean value accurate. The blue line indicates the calculated pdf and the red
dots stand for the actual computed eigenfrequencies.

4

CONCLUSION AND OUTLOOKThe combination of the finite element method, the infinite element method and spectral expansions is a promising tool to evaluate the eigenfrequencies of meta-atoms. That procedure is also applicable to bigger structures,
like arrays of meta-atoms or metamaterials. Hence, the Bragg type effects as well as the diffraction at theedges of the structure can be identified by the quality factor due to the eigenvalues being complex numbers. A
possible next step is the analysis of a two dimensional finite array of meta-atoms with random parameters, thatlead to a more complex system of equations, thus, more computational power. Another one is the identification1967

of the random parameters of 3D printed meta-atoms and the comparison between simulation and experimentconsidering the eigenfrequencies or the insertion loss.

Figure 4. First 4 resonant frequencies evaluated in the half-space formulation.

Figure 5. Probability density function of the first 4 resonant frequencies.

1968

REFERENCES[1] Peiró-Torres, M. P.; Redondo, J.; Bravo, J. M.; Sánchez Pérez, J. V. Open Noise Barriers Based on Sonic
Crystals, Advances in Noise Control in Transport Infrastructures, Transportation Research Procedia, 18,
2016, 392-398[2] Sánchez-Dehesa, J.; Garcia-Chocano, V. M.; Torrent, D.; Cervera, F.; Cabrera, S.; Simon, F. Noise control
by sonic crystal barriers made of recycled materials, The Journal of the Acoustical Society of America, 129(3), 2011, 1173-1183[3] Fredianelli, L.; Del Pizzo, A.; Licitra, G. Recent Developments in Sonic Crystals as Barriers for Road
Traffic Noise Mitigation, Environments, 6, 2019, 14[4] Chalmers, L; Elford, D. P.; Kusmartsev, F; Swallowe, G. M. Acoustic band gap formation in twodimensional locally resonant sonic crystals comprised of helmholtz resonators, International Journal of Modern Physics B - IJMPB, 23, 2009, 1-9[5] Elford, D. P.; Chalmers, L.; Kusmartsev, F.; Swallowe, G. M. Matryoshka locally resonant sonic crystal,
The Journal of the Acoustical Society of America, 130 (5), 2011, 2746-2755[6] Schot, S. Eighty years of Sommerfeld’s radiation condition, Historia Mathematica, 19 (4), 1992, 385-401[7] Astley, R. J. Infinite elements for wave problems: a review of current formulations and an assessment of
accuracy, International Journal for Numerical Methods in Engineering, 49, 2000, 951-976[8] Marburg, S.; Nolte, B. Computational Acoustics of Noise Propagation in Fluids - Finite and Boundary
Element Methods, Springer-Verlag, Berlin, Heidelberg (Germany), 2008
[9] Astley, R. J. Mapped spheroidal wave-envelope elements for unbounded wave problems, International Journal for Numerical Methods in Engineering, 41, 1998, 1235-1254[10] Astley, R. J.; Macaulay, G. J.; Coyette, J-P.; Cremers, L. Three-dimensional wave-envelope elements ofvariable order for acoustic radiation and scattering. Part I. Formulation in the frequency domain, The Journalof the Acoustical Society of America, 103, 1998, 49-63[11] Gerdes, K. The conjugated vs. the unconjugated infinite element method for the Helmholtz equation inexterior domains, Computer Methods in Applied Mechanics and Engineering, 152, 1998, 125-145[12] Stefanou, G. The stochastic finite element method: Past, present and future, Computer Methods in Applied
Mechanics and Engineering, 198, 2009, 1031-1051[13] Wiener, N. The homogeneous chaos, American Journal of Mathematics, 60, 1938, 897–93[14] Le Maitre, O.; Knio, O.M. Spectral Methods for Uncertainty Quantification, Springer-Verlag, Netherlands,
1, 2010[15] Xiu, D.; Karniadakis, G. The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations, SIAM
Journal on Scientific Computing, 24, 2002, 619-644[16] Hosder, S.; Walters, R.; Balch, M. Efficient sampling for non-intrusive polynomial chaos applications withmultiple uncertain input variables, Proceedings of the 9th AIAA Non- Deterministic Approaches Conference,
125, 2007, 1-16
[17] Moheit, L.; Marburg, S. Acoustic radiation modes and normal modes in exterior acoustic problems, The
Journal of the Acoustical Society of America, 141 (5), 2017, 3978-3978

1969

[18] Hammarling, S; Munro, C. J.; Tisseur, F., An Algorithm for the Complete Solution of Quadratic Eigenvalue
Problems, ACM Trans. Math. Softw., 39 (3), 2013, 18:1–18:19[19] Cheema, M.; Kirk, A., Accurate determination of the quality factor and tunneling distance of axisymmetricresonators for biosensing applications, Optics express, 21, 2013, 8724-35[20] Kao, C.-Y.; Santosa, F. Maximization of the quality factor of an optical resonator, Wave Motion, 45 (4),2008, 412-427

1970