Computation of dispersion diagrams for periodic porous materials modeled as equivalent fluids

The application of shift cell technique is presented and discussed for periodic porous media described with equivalent fluid models: as it can be found in literature, it consists in a reformulation of classical Floquet-Bloch (F-B) conditions, in which the phase shift of the boundary conditions, related to wave propagation, is integrated into the partial derivative operator. Consequently, the periodicity is included in the overall behavior of the structure, while continuity conditions are imposed at the edges of the unit cell. Its major advantage stands in allowing the introduction of a generic frequency dependence of porous material behavior, through the resolution a quadratic eigenvalue problem, providing an efficient way to compute the dispersion curves of a porous material modeled as an equivalent fluid. A validation and a computational cost comparison are performed between the shift cell technique and the classical F-B approach, pointing out that the first can provide, among its other advantages, a sensible computational time reduction for this kind of analyses. The derivation of the equivalent acoustic properties of the unit cell from its dispersion characteristics is also investigated. To this aim, group velocity matrix formulation and a branch-tracking algorithm are ∗Corresponding author Email address: dario.magliacano@univ-fcomte.fr, dario.magliacano@unina.it (Dario Magliacano ) Preprint submitted to Mechanical Systems and Signal Processing November 28, 2019 described. Some test cases are used for validating the proposed methodology.

described. Some test cases are used for validating the proposed methodology.
Keywords: vibroacoustics, porous material, shift cell, branch tracking, dispersion diagram, transmission loss List of symbols c 0 sound speed in the interstitial fluid C g group velocity E total energy E k kinetic energy E p potential energy I flow of energy j imaginary unit K bulk modulus k wave number in the material k 0 wave number in the interstitial fluid m mass p pressure p * conj(p) p 0 amplitude of the incident pressure r thickness of the unit cell S surface interested by incident pressure s side length v instantaneous local velocity v E energy transport speed x, y, z space variables Z 0 characteristic impedance of the interstitial fluid Z c characteristic impedance of the material Z s surface impedance of the material Γ domain boundary θ, φ angles of incidence ρ density of the material ρ 0 density of the interstitial fluid τ ∞ transmission coefficient Ω poro-elastic volume ω angular frequency

Introduction
The inclusion of vibroacoustic treatments at early stage of product development, through the use of porous media with periodic inclusions, is a powerful strategy for the achievement of lightweight sound packages and represents a convenient solution for manufacturing aspects [1]. 5 The main advantage of designing sound packages with periodic arrangements is that they can provide a combination of absorption effects, resonance effects and wave interference effects. These configurations can address different applications in transportation (aerospace, automotive, railway), energy and civil engineering sectors, where both weight and space, as well as vibroacous- 10 tic comfort, still remain as critical issues. Indeed, although porous materials are commonly used for vibroacoustic applications, they suffer from a lack of absorption at low frequencies compared to their efficiency at higher ones; this difficulty is usually overcome by multi-layering [2]. However, while reducing the impedance mismatch at the air-material interface, the efficiency of such 15 devices relies on the allowable thickness [1]. A more efficient way to enhance the low frequency performances of sound packages consists in embedding periodic inclusions in a porous layer [3,4] in order to create wave interferences or resonance effects that may play a positive role in the dynamics of the system. Therefore, numerical tools to 20 properly design sound packages are more and more useful. The classical approach, known as Floquet-Bloch (F-B) theory, is a branch of the theory of ordinary differential equations relating to the class of solutions to 1D periodic linear differential equations of the formẋ = A(t)x with A(t) a piecewise continuous periodic function with period T . Floquet's theorem, 25 due to Gaston Floquet [5], gives a canonical form for each fundamental matrix solution of this common linear system, through a coordinate change that transforms the periodic system to a traditional linear system with constant, real coefficients. In solid-state physics, the analogous result, extended to 3D systems, is known as Bloch's theorem [6]. In the literature dealing with 30 wave propagation problems in mechanical systems, the theory is referred to as Floquet-Bloch theory. In physical sciences and engineering, dispersion relations describe the effect of dispersion in a medium on the properties of a wave traveling within that medium, and therefore they offer a good perspective to explain the wave 35 field behavior inside bodies, relating the wavelength or wavenumber of a wave to its frequency. Dispersion may be caused either by geometric bound-ary conditions or by interaction of the waves with the transmitting medium and, in its presence, wave velocity is no longer uniquely defined, giving rise to the distinction of phase velocity and group velocity. For instance, the 40 Helmholtz equation is a known example of equation describing the spatial behavior: there, the physical periodic structure of the studied object translates into spatial periodicity of its coefficients. Therefore, the F-B theory can be applied to obtain the dispersive properties of different mechanical periodic systems, reducing the problem to the calculations performed in the so-called 45 unit cell under to certain specific boundary conditions derived from the F-B theory itself [7]. In order to develop efficient numerical techniques to handle the problem, the shift cell operator technique is herein considered. It allows the description of the propagation of all existing waves from the description of the unit 50 cell through the resolution of a quadratic eigenvalue problem. It belongs to the class of the k(ω) (wave number as a function of the angular frequency) methods, that allow computing dispersion curves for frequency-dependent problems, instead of using the classical ω(k) (angular frequency as a function of wave number) that leads to non-linear eigenvalue problems. Similar 55 techniques, which use a modified F-B approach in order to handle a k(ω) problem, can be found in literature [8,9,10]. The shift cell method consists in a reformulation of classical F-B conditions, in which the phase shift of the boundary conditions, together with the exponential amplitude decrease related to wave propagation, are integrated into the partial derivative op-60 erator; consequently, the periodicity is included in the overall behavior of the structure, while continuity conditions are imposed at the edges of the unit cell. Hence, this approach avoids condensation or non-linear eigenvalue solver, which are required by other ones. This technique has been successfully applied for describing the mechanical behavior of periodic structures embedding visco-elastic materials [11,12] or piezoelectric materials [13,14].
Here it is proposed an extension to equivalent fluid models of porous materials, which makes possible to overcome the limits of existing approaches in order to obtain a device whose frequency efficiency outperforms existing designs.

70
It should be pointed out that a homogeneous 3D unit cell with a 2D periodicity essentially represents an infinite layer with a given thickness; anyway in this paper, in order to keep a consistent nomenclature, this configuration is still addressed as "homogeneous case" or "homogeneous unit cell". Herein, the behavior of the porous material is described by Johnson-Champoux-Allard 75 4 (JCA) model [15,16], but one can identically use any other equivalent fluid model [17,18,19]. Then, starting from the complex wave numbers, obtained as an output from the quadratic eigenvalue problem, it is shown that it is possible to compute an equivalent transmission loss curve for an excitation at normal incidence; 80 the comparison with results obtained using classical methods, shows a very good agreement. The aim of this paper, therefore, is to introduce some enhancements to the state of the art of the shift cell technique applied to equivalent fluid models.

85
In this paper, a porous medium with a periodic arrangement of perfectly rigid inclusions is considered. The behavior of the foam is described by an equivalent fluid model in the frequency domain, i.e.: where p = p(x, ω) is the acoustic pressure, x = (x, y, z) is the coordinate vector, ω is the angular frequency, ρ = ρ(ω) is the equivalent fluid density 90 and K = K(ω) is the bulk modulus [2]. The periodicity is described by ρ(x − rn) − ρ(x) = 0 and K(x − rn) − K(x) = 0, ∀x ∈ Ω, where n is a vector of integers normal to the face considered, r = (r 1 , r 2 , r 3 ) is a matrix containing the three vectors defining the cell periodicity directions and lengths, and Ω is the domain of interest.

95
Eq. 1 applies everywhere except on the discontinuity surfaces, where appropriate boundary conditions apply. When finite densities and bulk moduli are used, these are the continuity conditions stated on pressure and normal velocity; when rigid inclusions are considered, the normal velocity on the inclusion surfaces vanishes.

100
The classical F-B approach, that here is recalled in its 1D formulation but can easily be generalized for 3D applications, provides a strategy to obtain a set of solutions of a linear ordinary equations system. Only the solution inside a period is needed, verifying that wheref is a solution of Eq. 1 and, according to the classical nomenclature, 105 β F = e k F L is called Floquet multiplier, while k F is the complex Floquet exponent. In addition, Floquet found that the solution at any point can be factorized in the following two terms: where f (x) is a periodic function, that represents the eigenvectors and carries the periodicity L of the coefficients of the problem [7]. Any solution of Eq.

110
1 can therefore be expressed in the form of Eq. 3. For the purpose of the shift cell technique development, considering Eq. 1 and applying the Bloch theorem, which generalizes Floquet's results to 3D systems, such as p(x, ω) = p(x)e jkx , where k, for a 3D application with real angles, is one can obtain p(x) being periodic, the Dirichlet boundary conditions imply continuity along the periodic directions. The results discussed in the following sections are obtained along the x-axis (φ = θ = 0 • ).

Weak formulation 120
The aim of this section is the development of the weak formulation of the problem, in order to obtain a matrix equation that fully describes what happens inside a periodic unit cell of equivalent fluid. A weak formulation of Eq. 5 consists in finding p such that ∀ p, which obeys to the periodic boundary conditions, one has: The solution approach follows a common weak formulation of a differential problem in a discrete coordinate scheme. After rewriting the second term through the use of an integration by parts, for which the considerations on classical weighted residual methods [20] are valid, and considering that Γ is the boundary domain, one obtains The boundary condition causes the integral on the boundary to vanish: Finally, one can discretize the weak formulation through the Finite Element Method: considering that ϕ is the eigenvector, the equation can be written in its matrix form with the following matrices: Here, M and K are respectively the symmetric mass and symmetric stiffness matrices, L is a skew-symmetric matrix and H is a symmetric matrix; all of them are complex and frequency-dependent.

Right and left eigenvalue problems
In this sub-section, the link between right and left eigenvectors is derived.

145
A left eigenvector of a matrix is the same as the right eigenvector of the same real transposed matrix. The formulation in Eq. 10 leads to the following right eigenvalue problem: where λ i = jk i is the i-th eigenvalue and ϕ r i denotes the right eigenvector associated to λ i . In this formulation, all matrices are frequency dependent.

150
The Eq. 11 can be rewritten as where I d is the identity matrix. Conversely, a left-eigenvector for the same eigenvalue satisfies In the resolution of the right eigenvalue problem, the i-th mode (i ∈ N + ) is defined by its λ i and its eigenvector ϕ r i . For each mode i, a mode −i is 160 associated with λ −i such that λ −i = −λ i and ϕ r −i = ϕ l i . Therefore, by solving the right eigenvalue problem, the left solution is found too [12].

Group velocity
For frequency-dependent systems, the estimation of the group velocity is not trivial [21]. In order to find its expression, Eq. 12 is now differentiated 165 with respect to ω: and multiplied by the left eigenvector such that: Considering that one obtains which gives the expression of the group slowness using λ i = jk i : 9 Finally, the complex group velocity is the inverse of the complex group slowness [21]: It should be noted that the frequency dependence of porous material models is generally known analytically, hence the computation of C g is fast since the derivatives of the matrices with respect to ω can be expressed analytically: with ∂ρ ∂ω to be derived from the specific equivalent fluid model. This is without doubts the most efficient way to perform the C g computation, but one could also choose to numerically estimate the matrix derivative, or even directly the ∆ω ∆k i values.

Classifying criteria to distinguish propagative and evanescent waves 180
For undamped systems, waves are classified according to their propagative (k purely real) or evanescent (k purely imaginary) behavior, but when dissipation occurs, such as it happens for a sample modeled through an equivalent fluid, all wave numbers are complex; consequently, there is no purely propagative solution and all waves are damped, with a decay rate that may 185 be used to classify the branches in two categories: those that are rapidly damped and those that are slowly damped in space. Hence, the latter could be classified as "propagative" ones. In general, the distinction between them is difficult and, thus, two classifying criteria are proposed.
I) The ratio between the real and the imaginary parts of every wave num- imag(k) . The physical meaning of C I is related to the fact that the real part of a wave number represents the propagative behavior, while its imaginary part is linked to the dissipation and therefore should be smaller than the real part in order to be able to consider a wave as propagative.

195
It should be pointed out that, since the real part of k is periodic while the imaginary one is not, in order to correctly apply this criterion, the real part of k must be turned into non-periodic, by mirroring it in correspondence to each period; in particular, starting with a = 2 and for each frequency f i with i > 0 of a specific dispersion branch, the 200 procedure is the following: II) The ratio between the real part of the energy transport speed, defined as v E = I E for undamped waves, and the real part of group velocity real(Cg) , where I is the flow of energy and E = E k + E p = 205 Ω 1 2 (ρv 2 + p 2 ρc 2 )dΩ is the total energy. Waves may be qualified as propagative when the energy is transported at a velocity which is at least close to the order of the group velocity.
Only the waves corresponding to C I > τ I and C II > τ II are considered as propagative ones. In practice, for the purpose of the following analysis, the 210 thresholds τ are chosen such as τ I = 1 and τ II = 0.7. These thresholds may be chosen differently according to the problem of interest [12]. Since there is no strict distinction between "propagative" and "evanescent" waves, an alternative would be to define an indicator of the "propagativeness" nature for each (ω, k) value of the dispersion diagram.

215
This is illustrated later in the results presented in Section 4.

Analysis of undamped case
In order to validate the shift cell technique implementation, in the studied configuration and for propagation along the x-axis, a first calculation is made 220 to compare shift cell results with those obtained using classical F-B periodic conditions, using (non-dissipative) air as material. Herein, all results are related to a 2D unit cell (top part of Figure 2), constituted by a square with side equal to 2 cm and with a 0.5 cm radius circular perfectly rigid inclusion located at the center of the unit cell, and to a 3D unit cell (bottom part of  It should be pointed out that, in this paper, the 3D cases are in fact 2D ones solved with 3D meshing, not exploiting the possibility, from the geometrical point of view, of doing a 2D meshing. This choice is motivated by the fact that a 3D mesh actually captures the behavior along an additional direction respect to the 2D one, allowing to perform analyses for every combination of 235 angles φ and θ. The comparison shows a perfect agreement between the results of the two methods ( Figure 2). In particular, one can observe that this arrangement exhibits a band gap between 6000 Hz and 10000 Hz for waves propagating along x direction. ing frequency of the Bragg band gap decreases when the radius is increased, and at the same time the width of the gap is increasing.
In the previous sub-section, only the real part of the wave number is shown. Now both real and imaginary parts are shown, the latter being actually positive but shown as negative in the plots due to axis consistence. If no damping 250 is included in the model, k is either purely real, the wave is then propagative, or purely imaginary, the wave being then evanescent (left column in Figure 3). On the contrary, instead of using the adiabatic value (142 kPa) for the bulk modulus of air, one can artificially add a frequency-constant imaginary part 255 to it (142 + j12 kPa is used here for illustration), so that one can simulate a band gap behavior in presence of dissipation. The new dispersion curves are shown in the right column of Figure 3. Indeed, a complex bulk modulus prevents the presence of ideal band gaps in dispersion curves; one can clearly see that the gap is opening but, because of the damping, k x is no longer 260 purely imaginary: around the band gap, the slow branch with undamped material becomes a fast wave when damping is added and allows rapid and damped energy transportation inside the band gap. So, the real part of k x (being mirrored and turned into non-periodic as explained in the 1 st classifying criterion described in Sub-section 2.4) is not equal to π r anymore, but remains low (compared to the imaginary one), which means that the wave will be strongly spatially attenuated. Also, a very fast branch, for which the imaginary part goes out of the plot bound, can be observed in the first two cases in the right column of Figure 3.

Comparison of computational cost 270
It is now performed a computational cost comparison between the shift cell technique and the classical F-B approach, pointing out that the first can provide, among its other advantages, a sensible computational time reduction for dispersion analyses. Figure 4 and Figure 5 show a comparison of the computational cost, in terms of time and as a function of the mesh size, between 275 the shift cell and the F-B techniques. In particular, both eigenproblems are solved using 100 frequency steps; the 2D unit cell is meshed using triangular elements while, for the 3D geometry, tetrahedral elements are used. Both 2D and 3D geometries correspond to those shown in the previous sections and these results, in terms of computational times, are related to an undamped 280 case. As a conclusion, for the case of interest, the calculation cost is always lower with the proposed approach than the one required by the classical Floquet-Bloch technique. The gain is increasing when the number of elements of the finite element 285 model is increasing, which makes the technique attractive. The lower cost is attributed to the management of the boundary conditions, which is much more simple in the proposed methodology, where continuity instead of F-B periodic conditions are needed.

Branch-tracking algorithm
The previous Section has highlighted that a dispersion diagram can have several branches, which are associated to different wave characteristics. In a dispersion diagram there is a set of points, forming branches, that one may wish to connect and follow according to the nature of each branch.

295
Some solutions are proposed in literature, such as a MAC sorting criterion [11], but these methods require to store many data at every iteration. Instead, the group velocity constitutes a relevant indicator in order to follow the branches from a point of calculation to the next one [12]. The proposed technique consists in comparing a single group velocity value at a specific  frequency C g i (f ) with C g (f + ∆f ): from the group velocity associated to a starting point, the routine compares the initial C g i with all the group velocities at the next frequency step f + ∆f and a minimization is made in order to identify the point at f + ∆f to which is associated the closest value of C g . Then, this point is defined as the new starting one and so on, step by 305 step, the branch is identified. In order to better appreciate the behavior of each branch in the frequency range of study, in the following plots, dispersion and C g curves are also colorized with a scale of colors that indicates the level of "propagativeness": the value 0 means that the wave at that specific frequency is totally spatially attenuated, while the value 1 represents a prop-310 erly propagative behavior. In particular, considering the criteria discussed in Sub-section 2.4, if all of them are satisfied then the propagativeness value is equal to 1, otherwise it is calculated as the product between the results of the two classifying ratios divided by the correspondent thresholds. It should be pointed out that, if a specific criterion is satisfied, its contribution to the 315 estimation of the level of propagativeness is always equal to 1, even if its related ratio is larger. In other words:

Results
In this Section, for a 3D melamine unit cell, whose geometric data are re-320 ported in detail in Sub-section 3.1, some results are shown in terms of evanescent -propagative dispersion and group velocity diagrams and branchtracked dispersion and group velocity diagrams. These curves are obtained for both a homogeneous and a heterogeneous (with inclusion) 3D unit cell, whose porous material is modeled with JCA model ([15, 16]). The analyses 325 are carried out in the frequency range 0 -17000 Hz. This range, indeed, is interesting for acoustic applications and assures that the wavelength is much larger than the pore size, which is a necessary condition in order to use equivalent fluid models. The size of the inclusion is also large compared to the typical characteristic length that may be observed on a representative  Table 1: Acoustical parameters of the tested porous material.
elementary volume describing the macroscopic behavior of the porous material [22]. It is well known that the parameters of the equivalent fluid models can have a strong impact on the performances of the acoustic device [23], hence they should be determined in a confident way. In the current case, they have been 335 experimentally determined and are reported in Table 1.
In Figures 6 and 7, the two plots on the top and the center of each figure represent two different ways to show the classification between evanescent and propagative waves: on the top, the plot is directly obtained from the application of the criteria presented in Sub-section 2.4, while on the center 340 the same distinction is shown using a color scale of "propagativeness" as described in Section 4.1. For each dispersion diagram, three eigenvectors are reported (Figure 8). Only the real parts of the eigenvectors, in terms of acoustic pressure field, are shown, the imaginary parts being null due to the fact that they are propor-345 tional to their correspondent real parts. They are all extracted at the frequency of 8500 Hz and along the direction that conventionally corresponds to θ = φ = 0 • in the 1 st Brillouin zone. Their branches are ordered as follows: at increasing frequencies, the 1 st branch is represented by the first real part that reaches the unitary value, the 2 nd one is the second that sat-350 isfies this condition and so on. It can be noticed that they respectively act along the x-, z-and y-axis directions. As already stated in the introduction, a homogeneous 3D unit cell with a 2D periodicity essentially represents an infinite layer with given thickness: in particular, continuity conditions are applied along x-axis and z-axis, while sound-hard wall boundary conditions 355 are used on the surfaces orthogonal to the y-axis. Therefore, in the homogeneous case, the so-called eigenvectors basically show the pressure acoustic field inside the medium in correspondence of the first three eigenfrequencies. For what concerns the propagative -evanescent plots, one can notice that the 1 st mode propagates at almost all frequencies, the 2 nd one appears to be propagative starting from middle frequencies, while the 3 rd and 4 th ones are relevant only at high frequencies. For all tested configurations, the branchtracking algorithm is able to correctly classify the solutions, even in the presence of band gaps, branch-crossing or branch-veering phenomena, as it can be appreciated from Figures 6 and 7.

Computation of transmission loss from dispersion diagrams
Dispersion diagrams presented above help designers to understand the nature of the waves that can propagate in a sound package and the way they are attenuated on the basis of an infinite periodic arrangement of the unit cell. In this section, it is shown how these results can be used to estimate 370 the transmission loss at normal incidence for an acoustic package composed by a finite arrangement of 5 cells. This, in a first approximation, allows comparing the dispersion relations and the acoustical characteristics of the equivalent finite medium. For more complex cases, advanced homogenization techniques may be used [24,25]. 375 For a plane wave configuration, the transmission loss is computed in three different ways. I) Transfer matrix method [2] (homogeneous case): with T 11 T 12 T 21 T 22 = cos(kr) j sin(kr)Z c j sin(kr) Zc cos(kr) .
II) Full FEM with 5 cells (case with inclusion): where Π incident and Π transmitted represent the incident and transmitted power, respectively. For this configuration, the calculation is made using an implementation of the plane wave forced response of the periodic cell accounting for fluid loading [26].    terest, there are no other propagative dispersion branches to which is associated a mode that acts along the plane wave direction in the TL analysis. For a plane wave that acts along an arbitrary direction, a more complex formulation is required. While k is a direct output of the dispersion relation, the equivalent characteristic impedance is com-395 puted as Z c = √ Kρ, where the density ρ is obtained from the JCA model and the bulk modulus is calculated as K = ρ( ω k ) 2 . The 3 rd way of computation actually consists in a homogenization and several papers that deal with the link between this kind of techniques and Bloch waves can be found in literature [27]. 400 Concerning the case with the inclusion, one can notice that an improvement of transmission loss properties, with respect to the homogeneous case, is shown at all frequencies, in particular in correspondence of a peak at a fre-quency around 7 kHz, where a gain of about 15 dB is observed, and at high frequencies.

405
Note that, for the sake of comparison with the related dispersion curves, only their 1 st branch, identified through the previously defined algorithm, is meaningful due to the fact that the corresponding mode is the only one, between those considered here (which are lowest order and therefore the least attenuated modes), that is actually excited during these transmission loss 410 simulations. Note also that the dispersion branch taken into account is actually propagative, according to the previously defined classifying criteria, in the whole frequency range considered. Indeed, the TL improvement peak exactly corresponds to the frequency range of the 1 st branch of dispersion curves in which 415 the wave is strongly spatially attenuated. This is definitely encouraging, for the purpose of deriving the equivalent acoustic properties of the unit cell from its dispersion characteristics.

Conclusions
In this paper, the formulation of the shift cell technique is presented to-420 gether with some applications aimed at evaluating the dispersive properties of some foams modeled as a periodic porous medium. The benefit of the k(ω) formulation, in the extraction of the eigenvalues, is shown and commented. Since the numerical formulation for the foams leads all complex wave numbers, one of the key elements of the present research activities is the definition 425 of adequate classifying criteria. As matter of fact, there is no purely propagative solution and all waves are damped, with a decay rate that may be used to classify the branches in two categories: those that are rapidly damped and those that are slowly damped in space. The proposed criteria allow to distinguish between the evanescent and propagative waves traveling in the 430 material and, consequently, to derive the group velocity. The proposed technique is validated through a comparison with the application of the classical Floquet-Bloch periodic conditions to a melamine unit cell. It is also evidenced a remarkable reduction of the computational costs associated with the application of this shift-cell technique. 435 Furthermore, thanks to a branch-tracking algorithm, it is possible to compute equivalent transmission loss curves, which show a very good agreement with those obtained with classical methods. The next steps of this research work will involve the Biot model [28] for 2D 24 and 3D geometries.