Applied Acoustics 157 (2020) 107026 Contents lists available at ScienceDirect Applied Acoustics journal homepage: www.elsevier.com/locate/apacoust Periodically arranged acoustic metamaterial in industrial applications: The need for uncertainty quantification J. Henneberg a,c,⇑, J.S. Gomez Nieto c, K. Sepahvand c, A. Gerlach a, H. Cebulla b, S. Marburg c a Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany Technical University of Chemnitz, Professorship for Textile Technologies, Reichenhainer Str. 31/33, 09126 Chemnitz, Germany c Technical University of Munich, Chair of Vibroacoustics of Vehicles and Machines, Boltzmann Str. 15, 85748 Garching, Germany b a r t i c l e i n f o Article history: Received 16 April 2019 Received in revised form 15 August 2019 Accepted 1 September 2019 Keywords: Tuned resonator Stop band materials Uncertainty quantification Generalized polynomial chaos Stochastic FEM a b s t r a c t Structures with periodically arranged resonators have attracted attention within vibroacoustics recently. These periodic structures attenuate wave propagation in a defined frequency range. Engineers can tune the resonators to produce structures with a desired band gap behavior. The majority of the established stop band materials is produced by additive manufacturing today. While promising scientific results have been achieved with this method, it is not yet fit for industrial applications. In order to develop new manufacturing approaches and to close the gap between science and industry, the influence of uncertainties of periodic structures has to be quantified. In this study, the influence of geometrical uncertainties on periodic structures is investigated. The authors develop a finite element model of an epoxy plate with beam resonators for this purpose. In a parameter study, the influence of some parameters is identified. Afterwards, the behavior of stop band material with uncertain input parameters is studied using spectral stochastic methods. Having predefined probability density functions, the generalized polynomial chaos expansion is used to propagate the uncertainty of these parameters. The stop band behavior is accordingly represented by the generalized polynomial chaos having unknown deterministic coefficients. A collocation-based stochastic simulation is used to estimate the coefficients employing the deterministic finite element model as a black-box. The results show the necessity of regarding uncertainties in periodic structures. The performed stochastic simulations are suitable to define manufacturing tolerances for the production of stop band material. Ó 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/). 1. Introduction During recent years, acoustic metamaterials received growing attention due to their possibility to gain acoustic properties which can not be found in nature. In particular, the possibility to attain structures with band gap behavior is useful in applications with vibroacoustic issues. A sonic crystal by locally resonant material is presented by [1] where rubber coated spheres of lead are embedded in an epoxy matrix to attain a three dimensional lattice. Thus, high sound transmission loss can be achieved at certain frequencies. Flexural waves in thin plates with spring-mass resonators are investigated in [2]. It is shown that wave propagation can be attenuated by locally attached spring-mass resonators. The behavior of different excitations is investigated with regard to ⇑ Corresponding author at: Robert Bosch GmbH, Robert-Bosch-Campus 1, 71272 Renningen, Germany. E-mail addresses: Johannes.Henneberg@de.bosch.com, johannes.henneberg@ tum.de (J. Henneberg). the vibration transmission. Similar structures have been considered with regard to the sound transmission loss in [3]. They demonstrated that the frequency band of increased sound transmission loss in a metamaterial-based plate can be broadened significantly by replacing a single resonator in each unit cell with multiple smaller appropriately damped resonators. The effect of periodic, tuned resonators with regard to vibroacoustics is investigated in the studies [4,5]. The stop band behavior of a sandwich plate with a stepped resonator has been explained in [6]. It has been shown that the sound transmission of the sandwich plate is significantly reduced within the stop band of the structure. The studies [7,8] investigate the attenuation of structural vibrations in a duct by applying locally resonant metamaterials. It is shown that the combination of different resonator configurations in one structure can lead to combined stop bands if they are appropriately designed. In [9], an approach using topology optimization for acoustic metamaterial is presented. In [10], it is shown that resonant stop band material is suitable to reduce mechanical crosscoupling in phased array transducers. Another approach realizing https://doi.org/10.1016/j.apacoust.2019.107026 0003-682X/Ó 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 2 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 acoustic metamaterial is using a very thin membrane having adjustable concentrated masses, known as membrane-type acoustic metamaterials [11,12]. Further approaches to realize structures with band gap behavior are presented in literature like [13–15]. Summarizing the results found in literature, stop band materials which are achieved with periodic structures of tuned resonators are a suitable solution for vibroacoustic applications. However, the application of these structures in industrial sectors is quite small. One reason for this is the missing evaluation of the influence of uncertainties. One of the few works is presented by [16] addressing seismic metamaterials with the focus of robust-to-uncertainties optimal design. The influence of point defects in a phononic crystal Timoshenko beam is investigated in [17]. In [18], the influence of point defects in periodic crystal thin plates is investigated. It is shown that the defect modes existing in the first band gap depend strongly on the size of the point defect and the filling fraction of the system. In [19], the wave finite element method is employed to investigate a structure with one-dimensional periodicity which is locally perturbed. The applied approach is based on the method presented in [20]. However, the impact of parameter uncertainty, e.g. manufacturing tolerances for structures with twodimensional periodicity has not been investigated. Up to now, no possibility to specify tolerances for the manufacturing of periodically arranged acoustic metamaterial is discussed in literature. An effective method to overcome such drawback is to use stochastic analysis in the design stage, as introduced in this paper. Uncertainty quantification in such material can be investigated using two techniques, i.e. sampling and non-sampling based methods. Monte Carlo simulation plays a major role in the former method in which large number of samples of uncertain parameters are required to realize the structure responses [21,19]. Though even the method is straightforward to use, it suffers the lack of fast convergence for large uncertainty and is ineffective in terms of computational time for a large number of parameter realizations. The latter, in contrast, shows reasonable accuracy and computational efficiency for representation of parameter uncertainty. The method is based on a radically different approach, namely constructing the functional dependence expressed in terms of spectral decomposition with unknown coefficients and orthogonal polynomial basis [22–27]. Among various available non-sampling techniques, the generalized Polynomial Chaos (gPC) expansion method has received considerable attention, particularly for non-Gaussian random variables [28,29]. The gPC is able to represent uncertainties in parameters as an expansion of orthogonal polynomials of standard random variables. Here, in this paper, this method is applied to propagate uncertain input parameters in acoustic metamaterials. The results show that uncertainty in the geometry of resonators leads to a randomly distributed band gap. Based on these results, manufacturing tolerances of the periodic structure can be tolerated to attain a band gap within a certain frequency band. This work is structured as follows: initially, we focus on the theoretical background of wave finite element method. In Section 3, an introduction to the spectral discretization of uncertain parameters is given. In an explanatory example, the application of the gPC is demonstrated. The finite element model is presented in Section 4. The band gap behavior and the results of the uncertainty quantification are given in Section 5. Finally, a conclusion is drawn. gation in such structures. A possible method to obtain the dispersion relations is an approach where the Bloch theorem is applied. The wave propagation can be described in a so called unit cell which is repeated infinitely [30]. A detailed elaboration of the theory, known as wave finite element method, can be found in [31– 33]. For freely propagating, time harmonic waves with angular frequency x a response variable w can be expressed by w ¼ Weiðxt kx x ky yÞ ð1Þ ; with W describing the wave mode through the thickness of the structure and kx ; ky being the components of the wavenumber in x and y direction. This can be obtained by employing finite element methods. Thus, the motion is described by a finite number of generalized displacements q. In case of one-dimensional periodicity, q ¼ ½q1 q2 ŠT contains the generalized displacements at nodes 1 and 2, cf. Fig. 1(a). According to Bloch theorem, the displacements are related to each other by [31] q2 ¼ kx q1 ; ð2Þ with kx ¼ e ikx
rx and lx ¼ kx rx ð3Þ where r x and lx represent the length of the periodic lattice and the propagation constant in the direction of periodicity, respectively. Thus, the generalized displacement vector can be described dependent on q1 KR ¼ ½I kx IŠT : with q ¼ KR q 1 ð4Þ The generalized displacements are related to generalized forces in the absence of damping by  K x2 M q ¼ f;  K x2 M KR q1 ¼ f:  ð5Þ with K; M; andf being the stiffness matrix, mass matrix, and force vector, respectively. Inserting Eq. (4) into Eq. (5) leads to  ð6Þ Assuming zero forces at the node 1 KL f ¼ 0   KL ¼ I k x 1 I with ð7Þ gives a reduced eigenvalue problem with  KL K x2 M KR q1 ¼ 0:  ð8Þ This eigenvalue problem can be solved for x which is a function of
the propagation constant lx . The solutions for x lx characterize the free wave propagation in the periodic structure. This approach can be applied easily to two-dimensional periodicity, cf. Fig. 1 (b). Thus, the vector of generalized displacements changes to q ¼ ½q1 q2 q3 q4 ŠT : ð9Þ Consequently, the reduction matrices are given by  T KR ¼ I kx I ky I kx ky I and h i 1 1 1 1 KL ¼ I kx I ky I kx ky I ; ð10Þ with 2. Dispersion in periodic, infinite structures kx ¼ e In this section, the employed method to calculate the dispersion relations in periodic structures applying the Bloch theorem is described. The dynamic behavior of periodic structures can be described by dispersion relations to characterize the wave propa- In equivalence to defined by ikx
rx ly ¼ ky ry and ky ¼ e iky
r y : ð11Þ lx , the propagation constant in y direction is ð12Þ 3 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Fig. 1. Generalized displacements q at various points in case of (a) one-dimensional periodicity (b) two-dimensional periodicity, based on [31]. As x   lx ; ly is a function of the propagation constants lx ; ly in case of two dimensional periodicity, the eigenvalue problem has to be solved for any combination of the propagation vectors in order to obtain a dispersion surface. Due to the periodicity of the structure, only the first Brillouin zone (BZ) with lx ; ly ¼ ½ p; pŠ needs to be considered to identify band gaps [34]. Dependent on the symmetry of the unit cell, the calculation of the band gap can be reduced to the contour of a part of the BZ, the so called irreducible Brillouin zone (IBZ). Thus, the computational effort to identify band gaps can be reduced. In case of asymmetric unit cells, special attention must be shown to the definition of the IBZ as shown in [35]. The relevance for the present study is further discussed in Section 5.2. 3. Spectral discretization of uncertain parameters 3.1. Theoretical background The spectral discretization methods are the key advantage for the efficient stochastic reduced basis representation of uncertain parameters in finite element modeling. This is because these methods provide a similar application of the deterministic Galerkin projection and collocation methods to reduce the order of complex systems. In this way, it is common to employ truncated expansions to discretize random input parameters and structural response. The generalized Polynomial Chaos (gPC) expansion permits such a discretization. The idea is that an uncertain quantity of interest can be approximated by summation of polynomials with random independent variables and deterministic coefficients. Thus, a new dimension is generated by the uncertainty which is sought as a set of orthogonal functions. The classical gPC expansion employs Hermite polynomials with multidimensional random variables as a trial basis for the probability space to represent uncertainty. Using this approach, the random parameter X is expanded as, see [23] for more details ð13Þ i¼0 where ai are deterministic unknown coefficient functions and Wi are a set of orthogonal polynomials in terms of the random variable n with the orthogonality relation of h i Wi ; Wj ¼ E Wi ; Wj ¼ E W2i dij ¼ h2i dij   ai ¼ hX ðnÞ; Wi ðnÞi D E Wi ðnÞ2 ð14Þ i ¼ 0; 1; 2; . . . ð15Þ where hX ðnÞ; Wi ðnÞi ¼ E½ X ðnÞ; Wi ðnފ ¼ Z X ðnÞWi ðnÞf ðnÞdn ð16Þ X here, f is the probability density function (PDF) corresponding to the random variable n and X denotes the random domain. For practical simulation, the gPC is truncated to a finite number of terms, hereafter denoted by N, X ðnÞ  In the following section, the applied method for the spectral discretization of uncertain parameters is introduced. In the first part, the theoretical background is discussed briefly. Afterwards, an explanatory example on how to apply the gPC is given. 1 X X ðnÞ ¼ ai Wi ðnÞ in which hi is the norm of the polynomial. The unknown deterministic coefficients can be determined using the Galerkin projection as N X ai Wi ðnÞ and i¼0 N¼ ðn þ pÞ! n!p! 1 ð17Þ in which, p is the maximum order of the expansion and n denotes the random dimensionality. While the first few coefficients of parameters with small uncertainty range seem to be enough for accurate gPC representation, higher order expansion is used for an accurate representation of parameters having large uncertainty. Details of numerical simulation using the gPC expansion are given in many works, cf. [28,23,36,37,24,38]. 3.2. Explanatory example In this section, a basic example is given for uncertainty quantification employing the approach of generalized polynomial chaos expansion. The examination for this example as well as for the later presented example with periodic structures follows the workflow shown in Fig. 2. The unknown coefficients of the gPC are calculated using the non-intrusive method, namely employing the stochastic collocation points. It allows the use of a deterministic black box model. Thus, it is suitable to perform the gPC to a wide range of models. For this basic example, the first eigenfrequency of a square plate under completely free boundary conditions is considered as objective of interest. The calculation is done with the classical plate theory. This theory is discussed in detail in [39,40]. The first eigenfrequency f 1 can be calculated by f1 ¼ k2 2pL2 3 Eh 12qð1 ! m2 Þ ; ð18Þ including the dimensionless frequency parameter k, the edge length L, the Young’s modulus E, the plate thickness h, the volume density q and the Poisson’s ratio m. Considering a square plate with completely free boundary conditions and Poisson’s ratio m ¼ 0:3; k2 denotes to 13:49 for the first eigenfrequency [40]. In this basic 4 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Fig. 2. The scheme illustrates the workflow to perfom the gPC employing the collocation method. Any arbitrary model can be used as a black box model. Note: In order to keep the overview in this scheme, the notation of the index of summation as well as the lower and upper bounds of summation are neglected. example, only one uncertain input parameter, the Young’s modulus E, is considered for the calculation of the eigenfrequency. The edge length L ¼ 100 mm, the plate thickness h ¼ 2 mm, and the volume density q ¼ 2:7 cmg 3 are deterministic input parameters. The uncertain parameter follows a Beta distribution with the shape factors a ¼ 4 and b ¼ 4 and lower and upper limit being 63 GPa and 77 GPa, respectively. The truncation order of the polynomial is chosen to be two. Corresponding to the Beta distribution, Jacobi polynomial provides optimal orthogonality. In the given example, the n2 54. The set of polynomials is W1;1 ¼ 1; W1;2 ¼ 4n, and W1;3 ¼ 45 4 collocation points CP are obtained by the roots of the third order Jacobi polynomial, resulting in CP 1 ¼ 0; CP 2 ¼ 0:522, and CP 3 ¼ 0:522. In the given example, the unknown deterministic coefficients are identified using the Galerkin projection resulting in a1;0 ¼ 7
104 ; a1;1 ¼ 1:357
103 , and a1;2 ¼ 0. As in this example only one uncertain parameter is considered, only one function X 1 is present. Employing the input polynomial and the collocation points, three realizations R1 ¼ 70000; R2 ¼ 73655, and R3 ¼ 66344 are obtained. Each of the realization denotes a single input of the Young’s modulus for the deterministic black box model. Hence, the eigenfrequency is calculated three times and the results are stored in a solution vector Sol. The set of collocation points is introduced in w, a set which is obtained by the truncation of n polynomials W1 . . . Wn . In this way, the matrix A is obtained. In this example, A results in 2 1:0000 6 A ¼ 4 1:0000 1:0000 0 2:09 2:09 Solving the equation Sol ¼ Ab 1:25 3 7 1:82 5: 1:82 ð19Þ ð20Þ leads to three coefficients b10 ¼ 93556:19; b11 ¼ 1170:01, and b12 ¼ 10:41. Together with the previous defined polynomials, these coefficients lead to a result function f. The number of result functions corresponds to the number of objective functions which are obtained from a single calculation of the black box model. In this example only the first eigenfrequency is obtained from the model. Hence, one result function f 1 ¼ 9:36
104 þ 4:68
103 n 1:17
102 n2 is obtained. This function can now be sampled with the random variable n which follows the previously defined distribution. Finally, the probability density function (PDF) as well as the statistical moments of the resulting function can be calculated. Thus, the influence of an uncertain input parameter, in this example the Young’s modulus, is quantified. Finally, the results obtained from the gPC are compared to results from a classical Monte Carlo simulation, cf. Table 1 and Fig. 3. For the sake of brevity, the theoretical background of the gPC is not fully discussed in this work. In [23], a comprehensive introduction as well as a tutorial to this topic is given. The example shows the possibility to realize uncertainty quantification with a small number of realizations employing the generalized polynomial chaos expansion. It should be denoted that the convergence gets slower when employing a higher number a uncorrelated uncertain input parameters. However, compared to classical Monte Carlo simulations, in most of the cases a faster convergence is obtained anyway. 4. Finite element model The spectral discretization method explained in the previous section allows us to investigate the impact of uncertainties in geometrical parameters of periodic structures on the stop band via a finite element model. Therefore, the test cases unit cell, 2  2 super cell (SC), 3  3 SC, and 4  4 SC are considered. The generic model represents a UC. By applying the Bloch theorem, this structure is assumed to be repeated infinitely in x and y direction, cf. Fig. 4. As the unit cell is repeated infinitely, the analysis with uncertain inputs in this unit cell can be interpreted as analysis of 5 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Table 1 Comparison of the results between a generalized polynomial chaos expansion (gPC) and a Monte Carlo simulation (MC) with different numbers of realizations (R). Mean value [Hz] Standard deviation [Hz] gPC R = 3 MC R = 100 MC R = 10000 MC R = 1 mio. 93556 1560 93677 1532 93588 1562 93556 1559 Fig. 3. Comparison of the results between a gPC and a Monte Carlo simulation (MC) with different numbers of realizations (R). Fig. 4. Finite element models used for the uncertainty quantification in periodic structures. From left to right: unit cell, 2  2 super cell and 3  3 super cell. different specimens wherein one specimen all periodic structures are identical. This may be useful to identify the influence of manufacturing tolerances which are caused by processes repeated after producing one specimen. However, the changes due to the uncertainty have to be identical for each unit cell within the periodic structure. Another important matter is the consideration of uncertainties within one specimen. As the unit cell is assumed to be repeated infinitely, this is not a suitable model to investigate this question. Therefore, a couple of unit cells are modeled together and is then used as a new unit cell. In adoption to the study [17], this cell is called super cell. In a m  m super cell, m is called the size of the super cell. The unit cell consists of a plate like structure with an added beam resonator. Its geometry is equal to the example presented in [41,42]. The resonator radius and length are 3.5 mm and 5 mm, respectively. The plate thickness is denoted to 1 mm. The finite element model is meshed with second-order hexagonal elements. In order to obtain a model with a small number of degrees of freedom and little mesh dependency of the results, a mesh convergence study is carried out, applying the guidelines given in [43]. Thus, the finite element model, cf. Fig. 4, shows less than 2% error for the band gap when changing the element edge length to halve. The simulations are carried out with the commercial software Comsol Multiphysics [44]. An epoxy material with Young’s modulus E ¼ 3500 MPa, Poisson’s ratio m ¼ 0:33, and mass density q ¼ 1:14 cmg 3 is considered. The material model is assumed to have linear elastic behavior in accordance to Hooke’s law. The above mentioned FE-models are considered then as blackbox to investigate the impact of parameter uncertainty on the band gap. To this end, the band gap is considered as random output approximated using the gPC expansion with unknown deterministic coefficients. The spectral stochastic method is employed, cf. [37], where a set of collocation points are generated in the random domain. The FE-models are then executed for realizations of samples at each collocation point. Once the realizations are known, a least-squares minimization procedure is used to calculate the 6 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 unknown deterministic coefficients. This leads to a closed form for the band gap depending on uncertainties in the random input parameters. The results are presented in the following section. 5. Numerical results The influence of uncertainty in input parameters on the band gap is investigated. The parameters are assumed as random variables having a predefined PDF. In this study, the uncertainties in parameters are assumed to be characterized by beta distributions owing to the fact that the parameters are defined within a bounded positive range. The boundary values of uncertain input parameters and shape factors of the employed beta distributions, a and b, are given in Table 2. By adjusting the boundary values and shape factors, the distributions can be easily adopted to a wide range of uncertain geometrical parameters arising in manufacturing processes. The width and the position of the band gap are considered as the random outputs of the structures. They are approximated by the gPC expansions having unknown deterministic coefficients and multi-dimensional orthogonal basis. The basis is constructed from the tensor dyadic product of the individual orthogonal polynomials employed as basis of uncertain input parameters, cf. [23] for details. 5.1. Unit cell In a first step, the deterministic solution of the unit cell model is calculated. Fig. 5 shows the Bravais lattice with the contour of the IBZ CXMC and the dispersion curves in a band diagram. The calculation along the contour of the IBZ is suitable due to the symmetry of the unit cell. Referring to the nomenclature of crystallography, the unit cell of this study represents a p4mm lattice, for details see [45,46]. The probability to obtain the full band gap on the conTable 2 Range of input uncertain parameters with associated distribution shape factors Betaða; bÞ. Uncertain parameter Plate thickness Resonator length Resonator radius Boundary values [mm] distribution factors ½0:9; 1:1Š ½4:5; 5:5Š ½3:15; 3:85Š Betað4; 4Þ Betað4; 4Þ Betað6; 6Þ tour of the IBZ is denoted to 99.46% for a p4mm lattice as shown by numerical test in [35]. The center position of the band gap and the associated width are denoted as Bp and Bw , respectively, as shown in Fig. 5. The calculation is easily done by Bp ¼ Bub þ Blb ; 2 ð21Þ with Blb and Bub being the lower and upper frequency boundary of the band gap, respectively. The width Bw is defined by Bw ¼ Bub Blb : ð22Þ Caused by the beam resonator, the periodic structure exhibits a resonant stop band. The band gap center position is obtained at Bp ¼ 15370 Hz and with a band gap width Bw ¼ 3240 Hz for the deterministic model. In a parameter study, the influence of different geometrical parameters to the band gap behavior of the structure is analyzed. Therefore, plate thickness, resonator length, and resonator radius are varied in a range from 90% to 110% with regard to the nominal values. The results are shown in Fig. 6. It can be seen that the band gap behavior is sensitive to geometrical changes of the investigated parameters. As a result, the band gap center position shows an almost linear dependency for all three parameters within the investigated ranges. The relation between plate thickness and resonator radius, and band gap center position are characterized by a positive gradient of nearly the same value. In contrast, the resonator length and band gap center position are related by a negative gradient. A similar behavior is found for plate thickness and resonator radius when estimating the band gap width. However, in this case the influence of uncertainty in the resonator radius has higher impact then of the plate thickness. A particular relation is obtained for the sensitivity to uncertainty in the resonator length. Especially when reducing the resonator length, the influence decreases. Furthermore, an extremum is obtained at a band gap width of 3300 Hz. The PDF of outputs for the unit cell are given in Fig. 7. It seems that the PDFs for the band gap center positions are identical for the three parameters considered as single uncertain inputs due to the fact that they have the same magnitude of the gradient. Concerning the band gap width, the lower influence of the plate thickness compared to the resonator radius results in a lower standard deviation. As in case of the parameter study, a particular behavior is obtained Fig. 5. Bravais lattice with the contour of the irreducible Brillouin zone CXMC and the dispersion curves of the unit cell with highlighted band gap center position Bp and band gap width Bw . J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 7 Fig. 6. Parameter study to identify the influence of plate thickness, resonator radius, and resonator length to band gap width and band gap center position. The value are investigated in a range from 90% to 110% from the nominal value of the parameter. in the PDF for the band width Bw when the resonator length is considered as an uncertain input parameter. The accumulation at a band gap width of 3300 Hz is related to the extremum which is found in the parameter study. In a secondary investigation, a unit cell with two uncertain input parameters, resonator length and radius, is considered. The resulting PDFs show similar mean values for band gap center position and band gap width compared to the investigations with only one uncertain input parameter as given in Table 3. However, taking two uncertain input parameters into account results in an increased standard deviation for both, band gap center position and band gap width. Furthermore, the PDFs show nearly symmetrical behavior. In order to evaluate the results obtained with the gPC, two cases (resonator length and resonator length + radius) have been calculated with a direct Monte Carlo simulation. Therefore, 10000 realization are considered to get a converged result. The results from gPC and Monte Carlo simulation show good accordance with respect to the mean value, standard deviation, and form of the PDF, cf. Fig. 7 and Table 3. As mentioned above, the unit cell is assumed to be repeated infinitely. Thus, the results can be interpreted as an uncertainty quantification in periodic structures where the uncertainty arises between different specimens but not within one specimen. Taking this into account, a scenario where the thickness varies between each specimen seem most likely. Assuming a molding manufacturing process, uncertainties are induced by the closing and opening process between every fabricated specimen. However, also resonator length and resonator radius could change between each specimen depending on the manufacturing process. In such a case, the given model of a unit cell is suitable to quantify the influence of uncertainties to the band gap behavior. 5.2. Super cell In contrast to unit cells, SCs offer the possibility to consider uncertainties arising within the periodic structure. Here, different independent input parameters are used to quantify the uncertainty in band gaps of resonators. The SC in this form does not show any symmetry. Thus, uncertainties which appear within one specimen are estimated. Applying the Bloch theorem to the SC, it is repeated infinitely as well. When calculating the band gap behavior in a SC with uncertain input parameters, special attention must be paid to the definition of the applied BZ. It is shown that the possibility to observe a full band gap on the contour of the IBZ strongly depends on the symmetry of the lattice [35]. Generally, a SC can be considered as a new unit cell because periodic boundary conditions are applied to its boundaries. Hence, the findings from [35] must be taken into account. By introducing uncertain input parameters to the geometry of the resonators, the SC does not show symmetry anymore. Hence, it belongs to the plane crystallographic group p1. In this case, the presence of a full, omni-directional band gap can only be verified by calculating the entire area of the BZ. This leads to a much higher computational effort to identify full band gaps within the SC. Furthermore, in a SC a larger number of uncertain input parameters is considered. Thus, more realizations are required to solve the gPC expansion. Finally, the computational time to solve the gPC expansion of a SC increases tremendously compared to the unit cell. In order to identify positions and values of the extrema within one dispersion surface, the full area of the BZ is calculated for a SC with the most asymmetrical realization used in the gPC, cf. Fig. 9. This realization is identified by the highest deviation of the center of mass compared the SC with nominal values of geometrical parameters, cf. Section 4. The three-dimensional plot of the dispersion surfaces of this SC is shown in Fig. 8. Fig. 8 (a) shows the first 16 dispersion surfaces of the structure. A detailed view is given in Fig. 8 (b) in order to provide more insight to the band gap behavior of the structure. Compared to the previous discussed unit cell, a larger number of dispersion surfaces are obtained in the equivalent frequency range. This is the result of a higher number of resonators which results in a higher number of modes when solving the eigenvalue problem, cf. Eq. (8). Thus, the dispersion surfaces within a given frequency range are obtained. Cutting along the contour of the BZ, cf. Fig. 9, would result in a two dimensional band diagram, similar to the 8 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Fig. 7. Probability density functions of band gap width and band gap position for various uncertain inputs parameters in a unit cell. Inputs: (a) plate thickness, (b) resonator radius, (c) resonator length, (d) radius and length of resonator, (e) resonator length with Monte Carlo simulation, (f) radius and length with Monte Carlo simulation. Table 3 Mean value and standard deviation from uncertainty quantification in a unit cell. The values are given for the band gap position Bp and the band gap width Bw . For comparison, the results of a Monte Carlo (MC) simulation are shown besides the results from the generalized polynomial chaos expansion. Uncertain input plate thickness resonator radius resonator length resonator length & radius resonator length MC resonator length & radius MC Bp Bw Bp Bw Bp Bw Bp Bw Bp Bw Bp Bw Mean value [Hz] Standard deviation [Hz] 15357 3234 15383 3233 15382 3197 15406 3195 15376 3190 15394 3191 450 61 390 104 424 118 589 164 425 116 579 158 one shown in Fig. 5. However, this would consist of 16 dispersion curves. The band gap is obtained between the 12th and 13th dispersion surface for the 2  2 SC. Thus, the band gap is defined by the maximum of dispersion surface 12 and minimum of dispersion surface 13. The minimum and maximum of each dispersion surface are indicated in Fig. 9. They are located on the contour of the first BZ CXMYNXC shown in Fig. 9. This is due to the fact that the deviation from the symmetry criteria is small. Similar results are obtained for the 3  3 SC. A detailed investigation of the results shows that all extrema can be found calculating only the points C; X; M; andY. Based on these results, the band gap calculation for the realization of the gPC expansion is performed only at 4 points of the contour of the BZ C; X; M; andY. Accordingly, the computational time is reduced significantly. It has to be pointed out that the calculation of the dispersion relations at certain points of the BZ is suitable for the investigated structure with the given distributions of the uncertain input parameters. For different structures and different 9 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Fig. 8. Dispersion surfaces of a 2  2 super cell with most asymmetric realization used for gPC expansion. (a) Shows the first 16 dispersion surfaces of the super cell. (b) Detailed view on the 12th and 13th dispersion surfaces. Between these dispersion surfaces a omni-directional band gap is present. The coloring corresponds to the frequency. Fig. 9. Dispersion surfaces of a 2  2 super cell with most asymmetric realization used for gPC expansion. Minima and maxima are indicated by . (a) 12th dispersion surface, (b) 13th dispersion surface. Contour of the BZ is indicated by distributions of the uncertain input parameters, it is necessary to prove whether the calculation is required for the full area of the BZ, only along the contour of the BZ, or only at certain point of the contour of the BZ. The results for uncertainty quantification in a unit cell, a 2  2, 3  3 and 4  4 SC are shown in Fig. 10. The resonator length is considered as independent uncertain input. Thus, 4, 9, and 16 uncorrelated resonator lengths with a beta distribution are considered in the 2  2, 3  3 and 4  4 SC, respectively. The mean values and standard deviations are given in Table 4. Analyzing the band gap center position, the deviation of the mean value is less than 1% comparing the different cell types. However, the standard deviation decreases significantly from 2.76% to 0.69% comparing the results from the unit cell and the 4  4 SC. The band gap width shows different behavior when taking more uncertain input parameters into account. Thus, the mean value of the band gap width is reduced slightly. The standard deviation is decreased from 3.69% in the unit cell to 1.60% in the 4  4 SC. Especially, the shape of the PDF changes when taking more uncertain input parameters into account. In order to characterize the obtained output distributions, the 3rd and 4rd central moments are considered, cf. Table 5. Employing the Pearson model and , respectively. of probability density function, it is possible to distinguish different types of distributions, cf. [29]. In the investigated test cases, the skewness of the PDF reduced when taking more uncertain input parameters into account. However, all PDFs are identified having a Pearson type 1 distribution, which is known as a generalized beta distribution. The mean value of the band gap center position shows convergence already for a 2  2 SC. Regarding the band gap width no clear convergence can be found. However, the analysis of the 3rd and 4rd statistical moment shows negligible differences between the 3  3 SC and 4  4 SC. Thus, the 3  3 SC provides an appropriate approximation of the influence of uncertainties in the investigated case. 5.3. Guidelines to quantify uncertainties in periodic structures Summarizing the results of this section, the authors provide a guideline to perform uncertainty quantification in periodic structures. While the gPC expansion can be applied straightforward in a unit cell, attention must be paid when using a SC approach especially in case of symmetric unit cells. In this case, a deviation from symmetry criteria is introduced by uncertain parameters. Thus, the 10 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 Fig. 10. Probability density functions of band gap width and band gap center position for a unit cell (UC), 2  2, 3  3, and 4  4 super cell (SC) with 1, 4, 9, and 16 independent, uncertain resonator lengths, respectively. Table 4 Mean value and standard deviation from uncertainty quantification in a unit cell and in super cells with resonator length being an uncertain input parameter. The values are given for the band gap position Bp and the band gap width Bw . Uncertain input UC 2  2 SC 3  3 SC 4  4 SC Bp Bw Bp Bw Bp Bw Bp Bw Mean value [Hz] Standard deviation [Hz] 15382 3197 15467 3154 15470 3087 15473 3053 424 118 211 70 147 57 107 49 cell approach of Section 5.2, a possible stop criterion for the band gap center position Bp could be  < 0:1% with    E Bpðm¼iÞ E Bpðm¼i   ¼ E Bpðm¼i 1Þ 1Þ  : ð23Þ Consequently, the stop criterion is fulfilled with the 3  3 super cell. However, any other stop criterion would be possible and can be adopted to the requirements of a specific task. In order to reduced the computational effort, new methods for the fast computation of band gaps, like presented in [47] or [48], could be implemented. 6. Conclusion Table 5 3rd and 4rd standardized moment from uncertainty quantification in a unit cell and in super cells with resonator length being an uncertain input parameter. The values are given for the band gap position Bp and the band gap width Bw . Uncertain input UC 2  2 SC 3  3 SC 4  4 SC Bp Bw Bp Bw Bp Bw Bp Bw 3rd central moment [–] 4rd central moment [–] 0.158 1.425 0.052 0.327 0.048 0.650 0.032 0.544 2.479 4.749 2.847 2.928 2.937 3.676 2.958 3.629 calculation of the band gap can not be reduced to the contour of the BZ as in a symmetric unit cell. The dispersion surfaces must be calculated for the first BZ. In order to reduce the computational effort, this should be done for the most asymmetrical realization required to perform the gPC expansion. In case the extrema to define the band gap are still located on the contour or at some characteristic points, the band gap calculation can be reduced for all realizations required. The workflow is given in Fig. 11. The size m of the super cell depends on the required approximation quality. Thus, some sort of stop criterion is defined. The size m is increased until the stop criterion is fulfilled. In case of the presented super A study on uncertainty quantification in periodic structures has been presented. The investigated structure consists of a plate like part with an added beam resonator. Thus, a band gap is achieved where no free wave propagation is possible. It has been characterized by the center position and the band gap width. In stochastic analysis, the influence of uncertainties in the geometrical parameters plate thickness, resonator radius and resonator length has been investigated. The spectral method based on the gPC is employed to represent uncertainty in parameters as well as in the center position and the width as outputs of finite element simulation. As a results, the presented approach offers a possibility to quantify the influence of uncertainties in periodic structures. The major contributions can be stated as follows: (i) the gPC expansion is adequately accurate to quantify the influence of uncertain geometrical input parameter in periodic structures, (ii) super cells consisting of multiple unit cells are suitable to approximate the influence of uncertainties arising within the periodic structure, (iii) band gap calculation in super cells consisting of symmetrical unit cells can be carried out at certain points of the BZ if deviation from symmetry criteria caused by introduction of uncertainties is low. The presence of band extrema on the BZ contour should be checked prior to the uncertainty quantification. For that, the dispersion surfaces within the full BZ area J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 11 Fig. 11. Guideline for uncertainty quantification in periodic structures for (a) unit cell approach to consider uncertainties between different specimens (b) super cell approach to approximate uncertainties within one specimen. are calculated for the realization with the highest deviation from symmetry criteria, and (iv) guidelines to apply uncertainty quantification in periodic structures are developed. The presented approach contributes to close the gap of uncertainty quantification in periodic structures. The obtained PDFs and statistical moments offer a feasible solution to define manufacturing tolerances for periodic structure in order to obtain a band gap behavior within certain requirements. This is important to transfer the promising results shown in scientific studies into industrial applications. Acknowledgement The authors gratefully acknowledge Robert Bosch GmbH Germany for the funding of this research. References [1] Liu Z, Zhang X, Mao Y, Zhu YY, Yang Z, Chan CT, Sheng P. Locally resonant sonic materials. Science 2000;289(5485):1734–6. [2] Xiao Y, Wen J, Wen X. Flexural wave band gaps in locally resonant thin plates with periodically attached spring-mass resonators. J Phys D: Appl Phys 2012;45(19): 195401. [3] Xiao Y, Wen J, Wen X. Sound transmission loss of metamaterial-based thin plates with multiple subwavelength arrays of attached resonators. J Sound Vib 2012;331(25):5408–23. [4] Claeys CC, Vergote K, Sas P, Desmet W. On the potential of tuned resonators to obtain low-frequency vibrational stop bands in periodic panels. J Sound Vib 2013;323:1418–36. [5] Claeys CC, Sas P, Desmet W. On the acoustic radiation efficiency of local resonance based stop band materials. J Sound Vib 2014;333:3203–13. [6] Song Y, Feng L, Wen J, Yu D, Wen X. Reduction of the sound transmission of a periodic sandwich plate using the stop band concept. Compos Struct 2015;128:428–36. [7] Claeys C, de Melo Filho NGR, Belle LV, Deckers E, Desmet W. Design and validation of metamaterials for multiple structural stop bands in waveguides. Extreme Mech Lett 2016:7–22. [8] Melo NF, Claeys C, Deckers E, Pluymers B, Desmet W. Dynamic metamaterials for structural stopband creation. SAE Int J Passeng Cars – Mech Syst 2016;9:1013–9. [9] Lu L, Yamamoto T, Otomori M, Yamada T, Izui K, Nishiwaki S. Topology optimization of an acoustic metamaterial with negative bulk modulus using local resonance. Finite Elem Anal Des 2013;72:1–12. [10] Henneberg J, Gerlach A, Storck H, Cebulla H, Marburg S. Reducing mechanical cross-coupling in phased array transducers using stop band material as backing. J Sound Vib 2018;424:352–64. [11] Ma G, Yang M, Xiao S, Yang Z, Sheng P. Acoustic metasurface with hybrid resonances. Nat Mater 2014;13(9):873–8. [12] Langfeldt F, Riecken J, Gleine W, von Estorff O. A membrane-type acoustic metamaterial with adjustable acoustic properties. J Sound Vib 2016;373: 1–18. [13] Cui ZY, Chen TN, Chen HL, Su YP. Experimental and calculated research on a large band gap constituting of tubes with periodic narrow slits. Appl Acoust 2009;70(8):1087–93. [14] Liu M, Xiang J, Zhong Y. The band gap and transmission characteristics investigation of local resonant quaternary phononic crystals with periodic coating. Appl Acoust 2015;100:10–7. [15] Lou J, He L, Yang J, Kitipornchai S, Wu H. Wave propagation in viscoelastic phononic crystal rods with internal resonators. Appl Acoust 2018;141:382–92. [16] Wagner P-R, Dertimanis VK, Chatzi EN, Beck JL. Robust-to-uncertainties optimal design of seismic metamaterials. J Eng Mech 2018;144(3):04017181. 12 J. Henneberg et al. / Applied Acoustics 157 (2020) 107026 [17] Chuang K-C, Zhang Z-Q, Wang H-X. Experimental study on slow flexural waves around the defect modes in a phononic crystal beam using fiber bragg gratings. Phys Lett A 2016;380(47):3963–9. [18] Yao Z-J, Yu G-L, Wang Y-S, Shi Z-F. Propagation of bending waves in phononic crystal thin plates with a point defect. Int J Solids Struct 2009;46(13): 2571–6. [19] Mencik J-M, Duhamel D. A wave finite element-based approach for the modeling of periodic structures with local perturbations. Finite Elem Anal Des 2016;121:40–51. [20] Mencik J-M, Duhamel D. A wave-based model reduction technique for the description of the dynamic behavior of periodic structures involving arbitraryshaped substructures and large-sized finite element models. Finite Elem Anal Des 2015;101:1–14. [21] Orris RM, Petyt M. Random response of periodic structures by a finite element technique. J Sound Vib 1975;43(1):1–8. [22] Ghanem RG, Spanos PD. Stochastic finite elements: a spectral approach. Mineola: Dover; 1991. [23] Sepahvand K, Marburg S, Hardtke H-J. Uncertainty quantification in stochastic systems using polynomial chaos expansion. Int J Appl Mech 2010;2 (2):305–53. [24] Sepahvand K, Marburg S. Random and stochastic structural acoustic analysis. In: Hambric SA, Sung SH, Nefske DJ, editors. Engineering vibroacoustic analysis: methods and applications, Ch. 10. New York: John Wiley & Sons; 2016. p. 305–38. [25] Xiu D, Karniadakis G. Modeling uncertainty in flow simulations via generalized polynomial chaos. J Comput Phys 2003;187(1):137–67. [26] Matthies H, Brenner C, Bucher C, Soares C. Uncertainties in probabilistic numerical analysis of structures and solids – stochastic finite elements. Struct Saf 1997;19(3):283–336. [27] Sepahvand K, Scheffler M, Marburg S. Uncertainty quantification in natural frequencies and radiated acoustic power of composite plates: analytical and experimental investigation. Appl Acoust 2015;87:23–9. [28] Xiu D, Karniadakis GE. The wiener–askey polynomial chaos for stochastic differential equations. SIAM J Sci Comput 2002;24(2):619–44. [29] Sepahvand K, Marburg S. On construction of uncertain material parameter using generalized polynomial chaos expansion from experimental data. Proc IUTAM 2013;6:4–17. [30] Bloch F. Über die Quantenmechanik der Elektronen in Kristallgittern. Z Phys 1929;52(7–8):555–600. [31] Mace BR, Manconi E. Modelling wave propagation in two-dimensional structures using finite element analysis. J Sound Vib 2008;318(4):884–902. [32] Ruzzene M, Scarpa FL. A general FEM technique to model wave propagation in cellular periodic structures. Proc. SPIE 5053, smart structures and materials 2003: active materials: behavior and mechanics, vol. 5053. p. 414–22. [33] Langley R. A note on the force boundary conditions for two-dimensional periodic structures with corner freedoms. J Sound Vib 1993;167(2):377–81. [34] Brillouin L. Wave propagation in periodic structures: electric filters and crystal lattices. Dover Publications; 1953. [35] Maurin F, Claeys C, Deckers E, Desmet W. Probability that a band-gap extremum is located on the irreducible brillouin-zone contour for the 17 different plane crystallographic lattices. Int J Solids Struct 2018;135:26–36. [36] Sepahvand K, Marburg S. Stochastic dynamic analysis of structures with spatially uncertain material parameters. Int. J. Struct. Stab. Dyn. 2014;14(8). 1440029(15). [37] Sepahvand K. Spectral stochastic finite element vibration analysis of fiberreinforced composites with random fiber orientation. Compos Struct 2016;145:119–28. [38] Sepahvand K. Stochastic finite element method for random harmonic analysis of composite plates with uncertain modal damping parameters. J Sound Vib 2017;400:1–12. [39] Leissa A. Vibration of plates. NASA SP, Scientific and Technical Information Division, National Aeronautics and Space Administration; 1969. [40] Blevins RD. Formulas for natural frequency and mode shape. Krieger Publishing Company; 1995. [41] Wu T-C, Wu T-T, Hsu J-C. Waveguiding and frequency selection of lamb waves in a plate with a periodic stubbed surface. Phys Rev B 2009;79: 104306. [42] Collet M, Ouisse M, Ruzzene M, Ichchou M. Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems. Int J Solids Struct 2011;48(20):2837–48. [43] Langer P, Maeder M, Guist C, Krause M, Marburg S. More than six elements per wavelength: the practical use of structural finite element models and their accuracy in comparison with experimental results. J Comput Acoust 2017;25 (4):1750025. [44] COMSOL Multiphysics Reference Manual. 5th ed.; 2017. [45] Kittel C. Introduction to solid state physics. Hoboken, NJ: J. Wiley; 2005. [46] Radaelli P. Symmetry in crystallography: understanding the international tables, IUCr texts on crystallography. OUP Oxford; 2011. [47] Boukadia RF, Droz C, Ichchou MN, Desmet W. A bloch wave reduction scheme for ultrafast band diagram and dynamic response computation in periodic structures. Finite Elem Anal Des 2018;148:1–12. [48] Casadei F, Rimoli J, Ruzzene M. Multiscale finite element analysis of wave propagation in periodic solids. Finite Elem Anal Des 2016;108:81–95.