Application of fractional order operators to the simulation of ducts with acoustic black hole terminations

Journal of Sound and Vibration 465 (2020) 115035

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Special Issue: Recent Advances in Acoustic Black Hole Research

Application of fractional order operators to the simulation of ducts with acoustic black hole terminations John P. Hollkamp, Fabio Semperlotti ∗ Department of Mechanical Engineering, Ray W. Herrick Laboratories, Purdue University, 177 S Russell St, West Lafayette, IN, 47907, USA

article info

abstract

Article history: Received 15 February 2019 Revised 14 October 2019 Accepted 18 October 2019 Available online 23 October 2019 Handling Editor: Victor V Krylov

We investigate the potential of fractional order operators to develop accurate and efficient models for acoustic black hole (ABH) terminations in ducts. While previous research has successfully developed methodologies to estimate the reflection coefficient of ABH terminations, the integration of ABH elements in larger structural components and assemblies often results in computationally intensive models. The size of the model is typically driven by spatial discretization requirements necessary to capture the power-law geometric variation required to achieve the ABH effect. In an effort to explore modeling approaches that could lead to effective simulations of structures with large numbers of embedded ABHs, we develop onedimensional fractional order models capable of emulating the absorbing behavior of an ABH termination in an acoustic duct. Two modeling approaches will be discussed based either on a fractional-order boundary condition or on a fractional-order acoustic domain of an ABH termination. In its most general form, our methodology synthesizes integro-differential models based on operators having complex-valued, fractional, and frequency-dependent order. The role and implications of these models for simulating structures with ABH elements are investigated and compared with the results from traditional modeling techniques in order to assess the viability of fractional order modeling for ABH systems. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Acoustic black hole Fractional calculus Homogenization Transfer matrix Complex order

1. Introduction The term acoustic black hole (ABH) indicates structural features embedded in either acoustic or structural waveguides that are capable of completely absorbing an incoming wave. The dynamic behavior of ABH elements stems from their characteristic geometry that follows a power law profile. As an acoustic (or elastic) wave enters the power law taper, both the phase and group velocities are progressively reduced while the amplitude of the particle displacement increases. In an ideal ABH (i.e. an ABH whose cross-sectional area vanishes at the end of the taper), the incoming wave can never reach the end of the taper; therefore, it is never reflected back. This intriguing behavior as an ideal wave absorber has garnered much attention particularly for applications to passive vibration and sound attenuation. We should also note, however, that in practical implementations the cross-section of the waveguide can never be reduced to zero. Therefore, some level of reflection should always be expected unless the ABH is combined with a damping mechanism. From a general perspective, acoustic black holes have been developed for both structural and acoustic applications. In structures, ABHs have been designed mostly to absorb flexural waves in beams and plates. Mironov [1] and Krylov et al. [2,3] pioneered the theoretical and design approaches for ABHs embedded in beams and plates. In this case, ABHs consist of either

∗ Corresponding author. E-mail address: [email protected] (F. Semperlotti).

https://doi.org/10.1016/j.jsv.2019.115035 0022-460X/© 2019 Elsevier Ltd. All rights reserved.

2

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

embedded tapers or slots having a power law thickness profile [1–8]. In acoustic waveguides, ABHs allow the absorption of pressure waves and are typically used to control pressure levels in ducts [9–12]. In this case, the power law decay of the duct radius as well as the variation of the wall admittance along the duct can be used to achieve the ABH behavior. This paper will specifically focus on the modeling of sound pressure wave propagation in air-filled ducts. A key objective of absorber designs in acoustic ducts is to achieve a low reflection coefficient across a broad frequency range. Before the consideration of ABH features, absorbers were typically designed by adding external damping. Walther [13] for example added a lossy medium with variable material properties to a duct to create acoustic absorbers whose impedance varied as a function of distance. This led to a low reflection coefficient across a wide frequency spectrum. Rather than relying on the addition of damping, an ABH utilizes the variation in the cross section geometry to produce the same broadband effect of a low reflection coefficient. Influential work in the modeling of acoustic black hole features in air-filled ducts was conducted by Mironov et al. [9] and, more recently, by Guasch et al. [10]. Both authors considered a duct having an ABH termination achieved by using a distribution of annular rings whose inner radii decreased following a power-law profile. In order to study this system, Mironov et al. [9] proposed an analytical approach leading to a closed form solution for the case of linear taper (m = 1), yielding both the pressure field in the duct and the reflection coefficient of the ABH termination. Later, Guasch et al. [10] extended the analysis of acoustic ducts to obtain the analytical solution and reflection coefficient for a quadratic (m = 2) ABH profile. Additionally, El Ouahabi et al. [11,12] built cylindrical acoustic ducts where the radius of the inserted rings varied linearly or quadraticly based on the work of Ref. [9] and measured the reflection coefficient to illustrate the black hole effect. Note that, unlike the case of structural elements, in acoustic waveguides the ABH behavior exists for taper exponents lower than m = 2 due to the combined effect of the variation of the wall admittance (for more details see Eq. (11) in Ref. [9]). For cases employing tapers with power law exponents m such that 1 < m < 2, analytical solutions are not available. In many applications involving ABH terminations (and, in some cases, also ABH embedded tapers), models are typically required to evaluate the dynamic response of the host system while the detailed dynamic response within the actual ABH region is often of secondary importance. Finite element techniques have been traditionally used to perform detailed numerical simulations of ABH systems [6,8,14]; however, the need to capture the power law thickness variation dictates minimum requirements on the model’s dimensionality which ultimately results in computationally intensive models. As an example, a thin plate with ABH terminations on the edges will require a full three-dimensional model instead of a more efficient Kirchhoff plate formulation. Computational efficiency becomes an even more critical requirement when considering systems with multiple embedded ABHs (e.g. structures with periodic ABH lattices [7]). The previous observation suggests that there is a need to identify methodologies capable of capturing the effect of the ABH on the dynamics of the host structure without requiring detailed modeling of the ABH itself. In an effort to explore methodologies for accurate and efficient dynamic simulations of structures with embedded ABH features, we consider the use of fractional calculus (FC). Fractional calculus is a branch of calculus involving integrals and derivatives of non-integer order. Although the field of FC (see Refs. [15–17] for detailed mathematical discussion) started in the seventeenth century, practical engineering applications have been quite limited. A fractional derivative of a function can be taken with respect to either time or space. A fractional time derivative is a natural tool to model memory-dependent systems [17–19]. Furthermore, the intrinsic damping nature of a fractional operator [20–23] allows time fractional derivatives to accurately represent dissipation in viscoelastic materials [24–26] and acoustic wave propagation in complex (mostly porous) media [27–29]. Thus, time fractional derivatives are inherently non-conservative operators. On the other hand, space fractional derivatives are indicative of non-local spatial effects as well as spatial attenuation in systems, such as fractal geometries [30–32], that potentially are still conservative. Our study explores the design and implementation of fractional order integro-differential models as an approach to capture the effective behavior of ABH features on the response of the host structure. The proposed approach will focus on space fractional derivatives as a means to capture the attenuation effect in ABH terminations. The underlying idea is that fractional order models could serve as effective (homogenized) models of ABH features. Hence, the fractional models capture the effect on the dynamics of the host structure without requiring the detailed simulation of the ABH geometry and dynamics. We present two fractional models, one using a fractional boundary condition (BC), and a second using a fractional domain. The methodology will rely on the value of the reflection coefficient calculated from the original ABH termination in order to determine the order of the equivalent fractional model. Ideally though, it would be preferable that the order of the fractional operator could be directly connected to the power law taper m. Although it might be theoretically possible to derive such a relationship, this study is intended to serve as a preliminary investigation into the application of fractional calculus to ABH based structures with the intent of stimulating the engineering community to combine these two concepts and explore the potentiality of FC. Nevertheless, it is certainly key to develop such a relationship in order to enable future developments of this methodology and to take full advantage of this approach in the context of practical applications. The remainder of this paper is organized as follows: first, we discuss the general case of an acoustic duct with an ABH termination and review the procedure to calculate the reflection coefficient. Next, we introduce the idea of fractional calculus to model ABH systems using a fractional boundary condition or a fractional domain. Once the fractional order formulations of the ABH are introduced, we use the models to analyze a specific duct configuration and present a numerical validation by means of finite difference solution of the governing equation.

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

3

Fig. 1. Depiction of an acoustic duct with an ABH termination. The drawing shows the top half of the axisymmetric system.

2. Calculation of the reflection coefficient due to an ABH termination In this section, we calculate the reflection coefficient of air-filled ducts with ABH terminations via the transfer matrix (TM) method. To accomplish this, we review and use the formulation developed in Refs. [9,10]. First, we provide the ABH duct model and its equation of motion as developed in Ref. [9]. From there, we implement the TM approach for the ABH duct, which was originally formed in Ref. [10]. Finally, this section will also present a parametric analysis of the effects of the discretization levels and of different taper profiles on the reflection coefficient. Note that the knowledge of the reflection coefficient represents a key step in order to synthesize equivalent fractional order models. 2.1. The acoustic black hole model and governing equation Consider the air-filled semi-infinite duct with an ABH termination as depicted in Fig. 1. The ABH termination was created using annular rings with power-law decreasing radii [9,10]. The rings’ thickness is assumed negligible. Starting from the linearized form of the continuity and momentum conservation equations, the pressure in the duct is described by Ref. [9] as

[

]

𝜕 2 p 𝜕 p 𝜕 (lnS(x)) 2Y (x) + + p k20 + jZ0 k0 =0 𝜕 x2 𝜕 x 𝜕 x h ( x)

(1)

where p is the pressure, S(x) is the cross sectional area of the duct, k0 = 𝜔∕c is the wavenumber, 𝜔 is the frequency, c is the speed of the acoustic wave, Z0 = 𝜌0 c is the characteristic impedance of air, 𝜌0 is the air density, h(x) is the radius of the cylindrical duct, Y(x) is the √ admittance of the wall, x is the axial coordinate with the origin located at the tip of the power-law taper (see Fig. 1), and j= Y (x) = (−j𝜔)

1 𝜌0 c2

−1. The wall admittance can be approximated by the continuous lumped admittance [9,10]

h20

− h ( x) 2 2h(x)

(2)

where h0 is the radius of the section of the duct having constant cross section. The radius of the cylindrical duct in the termination area is given by the power law equation h ( x) =

h0 |x − L|m Lm

(3)

where L is the axial length of the ABH and m is the power law exponent. Mironov et al. [9] presented the analytical solution of Eq. (1) and the reflection coefficient produced from the ABH termination for the case m = 1 while Guasch et al. [10] presented the analytical solution and reflection coefficient for the case m = 2. 2.2. Transfer matrix approach The analytical solution to Eq. (1) is not available for cases of non-integer taper coefficient m. For these remaining cases, a well accredited approach is based on the use of the transfer matrix method which provides an approximate solution for both the pressure field and the reflection coefficient. In the TM method, the ABH tapered section can be discretized in a series of shorter

4

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Fig. 2. Ideal versus discretized duct profile. In the TM approach, the ABH taper is approximated by a series of shorter ducts having constant cross sectional areas.

ducts with constant cross sections, as depicted in Fig. 2. This discretization approach is akin to using a rectangular Riemann sum approximation to estimate a definite integral. The transfer matrix for a duct of constant cross section is given by Ref. [33] as

[ ] Pi

Ui

𝜌 c ⎡ cos(k L ) j 0 sin(k0 LD )⎤ [ ] 0 D ⎢ ⎥ Po SD =⎢ S ⎥ D cos(k0 LD ) ⎥ Uo sin(k0 LD ) ⎢j ⎣ 𝜌0 c ⎦

(4)

where the subscripts i and o indicate the inlet and outlet sections of the duct, P is pressure and U is acoustic volumetric velocity, and LD and SD are the length and the cross sectional area of the duct. Recall that, in practical implementation, ABHs exhibit a residual thickness at their end (the ABH terminates at x = L − l in L −l Fig. 1 rather than x = L). By dividing the ABH domain into N ducts of equal length LD = , the cross sectional area of a generic N

section q is Sq = 𝜋 h2q where the radius hq is determined by taking the value of the radius from Eq. (3) at the midpoint of the discretized ducts (as exemplified in Fig. 2). Note that Eq. (4) does not account for the variation of the wall admittance along the ABH duct. This contribution could be accounted for by considering the impedance matrix [10]

[

]

1

0

(5)

Yqcav 1 where Yqcav = j

(

)

k0 1 𝜋 L h2 − (h2 + h2q + hq−1 hq ) . Z0 D 0 3 q−1

(6)

Incorporating the matrix of Eq. (5) in Eq. (4), the resulting transfer matrix for each duct subsection is given by

𝜌 c ] ⎡ cos(k0 LD ) j 0 sin(k0 LD )⎤ [ Sq ⎢ ⎥ 1 0 . Tq = ⎢ S ⎥ cav q 1 cos(k0 LD ) ⎥ Yq sin(k0 LD ) ⎢j ⎣ 𝜌0 c ⎦

(7)

The complete transfer matrix model describing the entire ABH duct can then be assembled as

[ ] Pi

Ui

[ ]

= Tt

Pl

(8)

Ul

where the subindices i and l represent quantities taken at the inlet and at the termination of the ABH. The complete transfer matrix Tt is

[

Tt =

T11t T12t T21t T22t

]

= T1 T2 T3 … .TN−1 TN .

(9)

For harmonic excitation, the pressure field in the duct (i.e. before the ABH termination) can be described by the plane wave solution p(x, t ) = Aej(𝜔t−k0 x) + Bej(𝜔t+k0 x)

(10)

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

5

Fig. 3. Plot of the reflection coefficient as a function of frequency for an ABH with taper coefficient m = 2. The plot is parameterized in terms of the spatial discretization parameter N.

where A and B are constants to be determined. The spatial profile of the pressure is P (x) = Ae−jk0 x + Bejk0 x . Dividing by A, we define R = B∕A where R is the reflection coefficient. Then, the pressure Pi at the inlet consists of an incident wave of unit amplitude and a reflected wave of amplitude R. That is Pi = 1 + R.

(11)

Likewise, using the linear Euler equation [9], the acoustic volumetric velocity at the inlet is Ui =

S

𝜌0 c

(1 − R)

(12)

where S is the cross sectional area at the inlet. Also, we assume that the edge termination of the ABH is completely rigid meaning that Ul = 0. Substituting this assumption along with Eqs. (11) and (12) into Eq. (8) yields 1 + R = T11t Pl S

𝜌0 c

(1 − R) = T21t Pl

(13) (14)

where Pl is the pressure at the rigid termination. These can accordingly be arranged into the matrix equation

[ ] ⎡ −1 T11t ⎤ R ⎡ 1 ⎤ ⎢ S ⎥ = ⎢ S ⎥. ⎢ ⎢ ⎥ T ⎥ P ⎣ 𝜌0 c 21t ⎦ l ⎣ 𝜌0 c ⎦

(15)

Thus, for a given value of 𝜔, one can solve for R using Eq. (15). 2.3. Effect of the duct discretization on the reflection coefficient Before proceeding to the fractional modeling of the ABH termination, we performed a numerical parametric study in order to evaluate the effect of the level of discretization (described by the number of divisions N) on the estimate of R based on Eq. (15). This is similar to the study in Ref. [10] where the effect of the number of rings on the estimate of R was analyzed. Consider an ABH taper with parameters h0 = 0.25 m, L = 0.5 m, l = 0.001 m, and m = 2. In order to contrast the reflection due to the edge truncation, we introduce losses by using a complex wave speed c = 340(1 + 0.05j). Fig. 3 plots the absolute value of the reflection coefficient R as a function of frequency. Four discretization cases (N = 10, 50, 100, and 500) are considered. As expected, Fig. 3 clearly shows that the value of R is very dependent on N. The smaller N is, the larger the oscillation amplitude of the reflection coefficient. As N increases, the trend converges to a smooth curve which coincides with the analytical expression derived in Ref. [10]. Given that the methodology presented below to retrieve an equivalent fractional order model of the ABH (see § 3, 3.1, and 3.2) relies on the target reflection coefficient obtained via the transfer matrix method, the accuracy of the target function R(𝜔) is critical. Hence, in the remainder of this paper, we consider a value N = 500 which provide a sufficiently smooth approximation of the tapered profile and, hence, of the reflection coefficient. From a practical point of view, the number and size of the rings will affect the black hole behavior of a real duct configuration. According to Ref. [10], “the high number of rings needed to recover the analytic ABHs may pose a severe limitation to practical

6

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Fig. 4. Plot of the reflection coefficient of the ABH as a function of frequency. The curves are parameterized based on different power law taper coefficients m. The range of m explored in this plot is consistent with the constraint m ≥ 1 provided by Mironov [9].

realizations of the ABH, which aim at a limited number of rings for manufacturing purposes.” It is likely that, for very thin rings, local resonances and structural modes will give rise to fluid-structure-interaction that will affect the performance of the ABH. Although Ref. [10] does investigate the effects of the size and thickness of the rings in his transfer matrix method, a series of experiments or very detailed fluid-structure-interaction simulations would be needed to further characterize the influence of the rings on the performance of the ABH. Nevertheless, it should be noted that the experimental work in Refs. [11,12] was able to demonstrate that a design of a finite number of rings could produce an acoustic black hole in a real configuration whose reflection coefficient was similar to the values obtained from the analytical formulas derived in Ref. [9]. 2.4. Effect of the power law exponent on the reflection coefficient In a similar fashion, we study the effect of different values of the taper exponent m on the reflection coefficient. This showcases that our discretized transfer matrix approach can still calculate the reflection coefficient for non-integer values of the taper exponent m. Consider the ABH where h0 = 0.25 m, L = 0.5 m, and l = 0.001 m, but m is non-integer. Fig. 4 plots the reflection coefficient parameterized with respect to the taper exponent m. Fig. 4 reveals some interesting trends. As m increases, the oscillation of the reflection coefficient is less drastic. In Ref. [10], the author notes that the quadratic taper usually produces a reflection coefficient smaller than the linear taper. However, this is not necessarily always the case. In fact, the curves in Fig. 4 corresponding to the non-integer values of m usually have a smaller reflection coefficient than the m = 2 case. Clearly, the optimal value of m in order to reduce the reflection coefficient depends on a variety of factors including the length of the taper, the location of the termination, and the frequency range of interest. 3. Fractional order models of an ABH termination It was previously highlighted that numerical modeling of systems with embedded ABHs can lead to very computationally intensive models due to the need to capture the thickness variation characteristic of an ABH. We address the need for a more effective modeling approach by using fractional calculus. More specifically, we explore two different approaches: 1) a boundary condition of fractional order, and 2) a fractional order element that replaces the ABH termination (see Fig. 5). In both approaches, space fractional derivatives will be used. Also, both models depend on knowledge of the reflection coefficient in order to determine the fractional parameters. As previously pointed out, this is a limitation of the current formulation but not of the approach in general. If an analytical relationship between the fractional parameters and the power law taper could be established, the fractional model could be used without the need for any other reference model. However, we highlight that, even in the current form, there might be situations where the fractional model presented here could be directly applicable. An example consists in the case where the reflection coefficient is experimentally determined [11,12]. In such case, the fractional model would serve as an effective homogenized model of the ABH duct. In the following sections, we detail the equivalent fractional model approach and develop the corresponding mathematics. Then, we apply both approaches to the simple single ABH example presented back in § 2.3 and compare the reflection coefficient obtained using the fractional model to that found via the TM method.

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

7

Fig. 5. Schematic representations of a) the traditional ABH model for a duct, b) the fractional order boundary condition model, c) the fractional order domain model used to represent an ABH.

3.1. The fractional boundary condition In the fractional BC approach, we replace the entire ABH tapered duct (Fig. 5a) with a fractional order BC that imitates the presence of the ABH termination (Fig. 5b). In this way, there is no need to solve for the actual response within the termination while the response in the main duct is still accurate because the correct equivalent acoustic impedance from the ABH is properly accounted for. The general form of a space fractional boundary is

𝜕p 𝜕𝛽 p | +c 𝛽| =0 𝜕t 𝜕 x |x=x0

(16)

where 𝛽 is the fractional order and c is the fractional wave speed; that is, a coefficient conceptually equivalent to the wave speed m𝛽

. For the specific analyses at hand, we set the amplitude of c to be equal in integer order systems but having dimensions of s to the speed of sound in air; further considerations on how to select this constant can be found in Ref. [34]. Note that Eq. (16) is similar in form to a typical one-dimensional absorbing boundary condition which is recovered exactly for 𝛽 = 1 [35,36]. The fractional derivative in Eq. (16) is taken as the left handed Riemann-Liouville derivative to mimic the effect the ABH would have on the dynamics of the main duct located to the left of the ABH. The left handed Riemann-Liouville derivative operator is defined as 1 dn d 𝛽 f ( x) = f (𝜉 )(x − 𝜉 )n−𝛽 −1 d𝜉 𝛽 dx Γ(n − 𝛽 ) dxn ∫a x

(17)

where Γ(x) is the Gamma function, n = ⌈𝛽 ⌉, and a is the lower terminal of the operator. 3.1.1. Identification of the equivalent fractional order In order for the fractional BC to exhibit a dynamic behavior equivalent to an ABH termination, we must determine the order 𝛽 in Eq. (16) such that the fractional BC will produce the same reflection coefficient as the tapered termination. Recall that the wave pressure in the main duct is given according to the ansatz p(x, t ) = ej(𝜔t−kx) + Rej(𝜔t+kx) 𝜔

(18)

where k = and R is the reflection coefficient. The reflection coefficient at a given frequency 𝜔 for a given ABH with power law c exponent m was determined in § 2.4 using the TM method. Substituting Eq. (18) into Eq. (16) (and taking the value a in Eq. (17)

8

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Fig. 6. a) The real part of the complex order 𝛽 as a function of frequency. b) The imaginary part of the complex order 𝛽 as a function of frequency. The curves are parameterized for different taper coefficients.

as −∞ to represent a semi-infinite structure) yields e−jkx0 + Rejkx0 − j

c

[

𝜔

] (−jk)𝛽 e−jkx0 + R(jk)𝛽 ejkx0 = 0

(19)

where x0 marks the location where the ABH taper begins. Since the coordinate system had been set up such that x0 = 0, Eq. (19) simplifies to 1+R−j

c

𝜔

[

] (−jk)𝛽 + R(jk)𝛽 = 0.

(20)

It can be seen from Eq. (20) that, when R = 0, the value of 𝛽 that satisfies the equation is 𝛽 = 1. Such a result is consistent with the general observation that when 𝛽 = 1, Eq. (16) reduces to the typical fully absorbing boundary condition (i.e. R = 0). Consider the reflection coefficients obtained in Fig. 4. To calculate the value of 𝛽 in Eq. (20) corresponding to R(𝜔), a numerical solver is used. Recall that the reflection coefficient is a complex quantity (only the magnitude is plotted in Fig. 4). The plot of 𝛽 (𝜔) obtained using the coefficients in Fig. 4 is given in Fig. 6. As expected, Fig. 6 shows that the tapers that exhibited larger oscillations of the reflection coefficient also exhibit larger oscillations in their corresponding fractional derivatives. Furthermore, the order 𝛽 necessary to match the reflection coefficient is frequency-dependent and is a complex number. The significance and meaning of the complex order is a challenging topic. In fact, the mathematics of complex fractional derivatives is still a relatively unexplored branch of fractional calculus. Even more challenging, applications involving complex fractional orders are rather sparse. Authors including Atanackovic et al. [37], Makris et al. [38], and Hollkamp et al. [34,39] have used complex fractional calculus to study various engineering systems. Despite the use of complex order fractional operators in these papers, the actual physical significance of a complex fractional order derivative is still not completely evident. Makris notes a link between complex-order derivatives to the modulation of both the amplitude and phase of harmonic components of a function. He shows that an “important difference between real-valued and complex-valued time derivatives is that phase modulation in the latter case is frequency dependent whereas in the former is not” [38]. Lastly, Re(𝛽 ) = 1 and Im(𝛽 ) = 0 is consistent with the fact that a zero reflection condition can be modeled by the usual integer-order absorbing boundary condition. Before moving on to verifying the validity and performance of the fractional boundary condition, we would like to provide some perspective about the general context for the use of fractional boundaries to represent ABH terminations. Certainly, for simple one-dimensional cases like the one under analysis, there are other approaches not involving fractional calculus that are capable of modeling partially reflected waves [33,40] and that could serve a purpose similar to the fractional order boundary. For example, one could replace the ABH taper with the integer order boundary condition (which we refer to as the integer-order partial reflection BC)

𝛾

𝜕p 𝜕p | +c | =0 𝜕t 𝜕 x | x =x 0

(21)

where 0 ≤ 𝛾 ≤ 1. The case 𝛾 = 1 represents zero reflection while 𝛾 = 0 indicates complete reflection. All the intermediate values of 𝛾 represent partial reflections. The main point of this section is to show how fractional order mathematics (here taken in its most simple form of a one dimensional boundary conditions) could be used for modeling ABH terminations and, more importantly, to present an example that serves as a gateway problem to stimulate future thinking for the application of this mathematical tool to acoustic and

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

9

vibration analysis. As will be more evident in § 3.2, fractional calculus presents many unique features that are particularly wellsuited to describe the behavior of acoustic problems in inhomogeneous media (e.g. porous media) that can find many interesting applications in ABH based systems. 3.1.2. Finite difference model of the acoustic duct and fractional BC In order to verify the validity and accurateness of the fractional BC approach, we perform a numerical simulation of a semiinfinite duct terminated by the proposed fractional order boundary. The governing equations for such a system can be solved numerically by finite differences (FD) and compared with the predictions from the integer order TM model in order to confirm that the structure of the reflected waves (hence the reflection coefficients) are indeed the same in both models. In the fractional boundary condition approach, the governing equation describing the pressure waves in the main duct is still the classical second order wave equation

𝜕2p 𝜕2p = c2 2 𝜕 t2 𝜕x

(22)

which, by using a second-order centered scheme, can be written in an equivalent finite difference form pbn+1 − 2pbn + pbn−1

Δt 2

= c2

pbn+1 − 2pbn + pbn−1

(23)

Δx2

where n = 0,1,2, …,N is a spatial index, b is the time index, and Δx and Δt describe the spatial and temporal discretization, respectively. If the initial conditions are known (corresponding to b = 0,1), Eq. (23) can be solved for pbn+1 and then subsequently be placed in a computational loop. Note that the explicit FD scheme in Eq. (23) is stable if cΔt < 1. Δx

(24)

In order to represent the fractional boundary condition in the computational simulation, we developed a unique discretization of the fractional BC based on a FD formula known as Mur’s BC [41]. In a typical one-dimensional wave problem, the Mur BC can be implemented to model the traditional fully absorbing boundary condition. We modified the Mur BC to account for the space fractional derivative in Eq. (16). We indicate this new form of the Mur BC as the fractional Mur BC. In order to formulate the fractional form of the Mur BC, we must first introduce the FD definition of the Riemann-Liouville fractional derivative, better known as the Grunwald-Letnikov definition of the fractional derivative. The Grunwald-Letnikov derivative can be defined recursively as [42] v 1 ∑ d 𝛽 f ( x) = g f (x ) dx𝛽 (Δx)𝛽 q=0 q q

(25)

g0 = 1,

(26)

where

( 1−

gq =

𝛽 +1

)

q

gq−1 .

(27)

Using the Grunwald-Letnikov definition, the fractional Mur boundary condition is (See Appendix A for details)

( pbN+1

=

pbN

1−r r+1

)

(

+

1 r+1

)

(

(pbN−1

−

pbN+−11 )

−

r r+1

) (∑ N q=1

gq pbN+−1q

+

N ∑

) gq pbN−q

(28)

q=1

where the subindex N represents the last spatial node and r=c

Δt . (Δx)𝛽

(29)

3.1.3. Numerical validation The discretized model developed in § 3.1.2 can be used to obtain the response of the duct with a fractional BC and to compare the response with the traditional model. More specifically, we consider a situation where a wave packet (i.e. a Hanningwindowed wave at a given frequency) travels in the main duct towards the fractional BC. After encountering the fractional BC, the wave is partially reflected back. If the fractional BC serves as a representative model of the ABH termination, then the reflected waves should be comparable. Note that, according to the current strategy, this simulation in the time domain will only be valid for an impinging harmonic wave possessing the frequency corresponding to that of the fractional order used. However, the wave equation and the fractional

10

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Fig. 7. Numerical time-domain simulations performed by solving the finite difference model for the fractional BC (§ 3.1.2). The results show the response in the duct due to an incident wave packet at a selected frequency of Ω = 2988 rad/s. a) Real and imaginary parts of the incident pressure wave. b) Real and imaginary parts of the reflected pressure wave field after the entire waveform has encountered the fractional Mur boundary.

boundary condition are linear, so it is anticipated that linear superposition will still hold. The current study is intended to illustrate the basic concept although further refinements are certainly needed to yield a more general approach. The modeled duct has the same numerical parameters given in § 2.3. The reflection coefficient and the fractional derivative associated with this ABH were shown by the black, solid curves in Figs. 4 and 6, respectively. For this numerical analysis we selected a frequency Ω = 2988 rad/s. While, in principle, any frequency could be used in this numerical test, the trend of both the reflection coefficient and the fractional order were particularly smooth at the selected frequency. After considering this specific frequency, we will consider the response across a wider frequency range. From Figs. 4 and 6, the frequency Ω = 2988 rad/s corresponds to a reflection coefficient of |R| = 0.1129 and a fractional order of 𝛽 = 1.074–0.05477j. These values are directly used in the discretized model. Other parameters for the discretization in both time and space were selected as Δx = 0.01 m and Δt = 1E-5 s. The incident wave packet was generated by windowing the inhomogeneous boundary condition p(0, t) = exp(jΩt) on the left end of the duct. This harmonic excitation on the initial spatial node persisted for a time of t = 20.5T where T is the duration of a period (T = 2𝜋 ∕Ω). A plot of the incident wave field before it encounters the fractional BC is given in Fig. 7a. As the incident wave encounters the fractional BC at the right end of the domain, it is reflected back as shown in Fig. 7b. Notice that the reflected pressure is a complex quantity just like the incident pressure field. If the incident wave had been only a cosine function, the real parts of the plots in Fig. 7 would correspond to the actual physical responses of the wave in the duct. Conversely, if the incident wave was given as a sine function, the imaginary parts of the plots in Fig. 7 would correspond to the actual physical responses. To measure the reflection coefficient with an incident wave represented using a complex exponential function, we calculate the ratio of peak values of the absolute magnitude of the incident and reflected wave packets from Fig. 7. Doing so yields a reflection coefficient of R = 0.1119. Recalling that the reflection coefficient using the traditional TM model shown in Fig. 4 was R = 0.1129 at Ω = 2988 rad/s, the percent error between the transfer matrix method and the fractional BC model is only 0.9%. However, due to the finite duration of the signal, the spectrum of the incident and reflected signals will not be single frequency. Thus, another potentially more accurate way to calculate the reflection coefficient is to analyze the spectrum of either the real or imaginary parts of the incident and reflected waves by taking the fast Fourier transforms (FFTs) of the signals and measuring the ratio of the amplitudes of the FFTs at the wavenumber of interest. This approach yielded R = 0.112. It is not surprising that this value is almost nearly equal to the value determined via the ratio of the peak values approach since the Hanning window of the signal was chosen so that contributions of other nearby frequency components were minimized. These results confirm the validity of our procedure as well as the satisfactory performance of the fractional order boundary as a means to simulate the effect of the ABH termination on the response of the main duct at the selected frequency of Ω = 2988 rad/s. For further validation, we extended the analysis to a wider range of frequencies. The results are shown in Fig. 8 by plotting the reflection coefficients from the TM method and the FD time domain simulation. For the FD simulation curve in Fig. 8a, the reflection coefficient was calculated from the ratio of amplitude of the FFT component at the corresponding wavenumber of interest. The ratio between the real part of the reflected and the incident wave signal was used. Note that taking the ratio of the FFTs of the imaginary part of the signals would yield the same reflection coefficient curve. The curves in Fig. 8a nearly lie on top of each other. The percent difference between the curves is plotted in Fig. 8b. The error is always less than 1%, further proving that the fractional BC is an accurate and appropriate means to model the reflection from the acoustic black hole. A supplementary analysis available in Appendix B obtained the results from a FD simulation which used the integer-order partial reflection BC (Eq. (21)) in order to compare to the results from the FD simulation using the fractional BC (Eq. (16)).

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

11

Fig. 8. a) Plot of the reflection coefficients from the TM method and the FD simulation over a frequency range from 450 to 3000 rad/s. The plot was cut at 450 rad/s due to large computational expenses of the FD simulation for very low frequencies. b) Plot of the percent difference error between the curves of the reflection coefficients. Sources of the error include the inherent error of the finite difference method and a small error in the calculation of a complex-valued fractional derivative.

3.2. Fractional domain model and fractional transfer matrix The second approach that we explore in this study consists in replacing the tapered ABH termination with a finite domain of fractional order (that is a domain described by a space fractional wave equation). This fractional domain has the same length as the ABH taper but a constant cross sectional area with a value equal to the area at the inlet of the ABH, as shown in Fig. 5c. In this approach, the spatial dependence of the parameters (namely the wave speed) due to the cross sectional variation in the ABH section is replaced by a domain possessing a constant cross section and a space fractional derivative. In other terms, the inhomogeneities due to the variable cross section are lumped into an equivalent space fractional order wave model giving rise to an equivalent homogenized model of the ABH. As for the fractional boundary case, the order of the fractional domain must be selected to exhibit the same reflection and transmission properties of the actual ABH. Contrarily to the previous case of a fractional boundary, the fractional domain provides a physical domain in which propagation can take place. However, we are not concerned in matching the dynamics of the wave inside the fractional domain to that of the ABH; instead, we want to achieve an equivalent response within the duct. We note that, even if interested in matching the dynamic response inside the ABH, a fixed fractional order derivative could not capture exactly the wave behavior within the ABH termination because the order would change in space as a function of the local thickness. In order to simulate the response of a duct terminated by an ABH-like fractional domain, we develop a fractional order transfer matrix model and follow a similar simulation approach as that used in § 2.2 for an integer order system. 3.2.1. Fractional order transfer matrix The governing equation of the fractional domain is the space-fractional wave equation [23,34] c2

𝜕𝛼 p 𝜕2p = 2 𝜕 x𝛼 𝜕t

(30) (m)𝛼 ∕2

. The solution to this where 𝛼 is the order of the fractional derivative and c is the fractional wave speed with dimensions s partial differential equation can be expressed using exponential functions if the lower bound of the integral in the fractional derivative (the value a in Eq. (17)) is −∞ [34]. Thus, the solution of the fractional differential Eq. (30) is ̂

̂

p(x, t ) = Aej(𝜔t−kx) + Bej(𝜔t+kx)

(31)

where ̂ k is a complex-valued wavenumber. This solution corresponds to forward and backward propagating waves, respectively. The relationship between ̂ k and the fractional order 𝛼 is obtained from the dispersion equation of Eq. (30) [23,34] as

(−ĵ k) 𝛼 = −

𝜔2 c2

.

(32)

The procedure to derive the transfer matrix for the fractional domain is very similar to that discussed above for an integer order model since both domain solutions are based on exponential kernels. As such, the fractional transfer matrix is

[

P1

U1

]

𝜌 c [ ] ⎡ cos(̂ kL) j 0 sin(̂ kL)⎤ P 2 S ⎢ ⎥ = ⎢j S sin(̂ kL) cos(̂ kL) ⎥ U2 ⎣ 𝜌0 c ⎦

(33)

12

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Fig. 9. Plot of the complex, fractional order 𝛼 of the fractional domain that replaces the ABH. Using a fractional domain with this order produces the same reflection coefficient as the ABH taper.

where P1 and U1 are the pressure and acoustic volumetric velocity of the left end of the fractional domain, P2 and U2 are the pressure and acoustic volumetric velocity of the right end of the fractional domain, L is the entire length of the fractional domain, k is a complex wavenumber. S is the constant cross-sectional area, and ̂

3.2.2. Identification of the equivalent fractional order In order to obtain the order 𝛼 of the ABH-like fractional domain, we must guarantee that the equivalent dynamic behavior of the duct is unchanged when using either a classical ABH termination or a fractional domain model. Equivalently, the two should exhibit the same reflection coefficient. k. Since the end of the ABH at x = L − l was assumed In order to determine the order 𝛼 , we need to determine the value of ̂ to be rigidly terminated U2 = 0.

(34)

Substituting this into Eq. (33), P1 and U1 can be written as a function of P2 . Eliminating P2 from these equations gives the P input impedance 1 of the fractional domain as U1

𝜌 c P1 = −j 0 cot(̂ kL). U1 S

(35)

If the starting point of the fractional domain is kept at x = 0, then P1 = 1 + R

(36)

and U1 =

S

𝜌0 c

(1 − R).

(37)

Substituting these two equations into Eq. (35) yields 1+R = −j cot(̂ kL). 1−R

(38)

The value of the reflection coefficient R of the ABH was obtained using the transfer matrix method in § 2. Substituting the value

of R into Eq. (38), a numerical solver can obtain ̂ k. Once ̂ k is obtained, Eq. (32) is used to find 𝛼 . Eq. (32) rearranged to solve for 𝛼 is

𝛼=

2 ln(𝜔∕c) − j𝜋 ln(−ĵ k)

.

(39)

This methodology can be tested by leveraging the same numerical example (§ 2.3) used in previous sections. Eqs. (38) and (39) can be solved to calculate the fractional order 𝛼 (Fig. 9) of the fractional domain at the selected frequency.

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

13

Fig. 9 reveals, once again, that the methodology produces fractional differential operators of complex and frequencydependent order. We observe that 1 < Re(𝛼 ) < 2 while Im(𝛼 ) is positive (except for very low frequencies). The real part of 𝛼 is responsible for the attenuating effect characteristic of the fractional domain [20–23]. At first, this may seem to contradict the expected behavior inside the ABH termination, since the wave speed should approach zero while the amplitude of the pressure should be amplified. Recall, however, that the model is not attempting to match the dynamic behavior within the ABH but instead, it focuses on the response of the main duct. To produce a reflected wave with the same reduced amplitude resulting from a classical ABH, the fractional domain itself must attenuate the amplitude of the pressure wave. 4. Conclusions This paper presented the results of an exploratory study focusing on the use of fractional order operators to simulate systems with embedded acoustic black holes. In particular, we studied the possibility to develop equivalent space fractional order models of ABH terminations in air-filled acoustic ducts. These models are intended to replace the ABH termination with an effective model capable of inducing an equivalent dynamic behavior of the pressure field in the main duct without requiring the full solution of the pressure field in the termination. The fractional order model can be interpreted as a homogenized model of the ABH termination. Two modeling approaches were presented. The first was based on a fractional order boundary condition, while the second on a finite fractional order domain. In both cases, the key factor was to determine the order of the fractional operator capable of producing a dynamic behavior in the duct equivalent to that produced by the real ABH termination. In other terms, the proposed approach does not focus on the determination of the detailed pressure field in the ABH termination, but instead focuses on replacing it with an equivalent model having comparable dynamic behavior and wave reflection characteristics. The method presented in this study uses the value of the reflection coefficient calculated from the original ABH termination in order to determine the order of the equivalent fractional model. While this approach currently restricts the applicability of the methodology (since the value of the reflection coefficient must be known a priori), this study shows the feasibility and the potential of such an approach. In summary, this work explored the application of fractional calculus to the simulation of one-dimensional acoustic systems with acoustic black hole terminations. The rationale behind such an approach is to limit the computational burden associated with the detailed modeling of the ABH geometry by developing effective (homogenized) models of an ABH capable of producing an equivalent dynamic behavior. Although this study concentrated on the development of one-dimensional models for air-filled acoustic duct, we anticipate that similar considerations and methodologies could be extended to more complex systems, such as one or two-dimensional structures with periodic distributions of acoustic black holes. Acknowledgments This work was supported by the United States National Defense Science and Engineering Graduate Fellowship, and by the National Science Foundation under MOMS #1761423. Appendix A. Derivation of the Fractional Mur Boundary Condition The fractional BC is rewritten as

𝜕𝛽 p 1 𝜕p =− . 𝜕 x𝛽 c 𝜕t

(A.1)

This partial differential equation is applied at the right end of the model. Following a derivation similar to that given in Ref. [41], we need to ensure that the spatial and time derivative approximations are evaluated at the same point. We do this by averaging the finite differences evaluated at the spatial nodes N and N − 1 and the temporal points b and b + 1. Thus, the fractional Mur boundary condition will be evaluated at the point (xN − Δx∕2, tb + Δt∕2). Note that while the Grunwald-Letnikov definition for the spatial fractional derivative considers all the spatial nodes in the domain, the values of the Grunwald weights are such that the evaluation of the spatial fractional derivative occurs at xN − Δx∕2. Using this procedure yields

[

1 2

1 Δx𝛽

(N ∑ q=0

gq pbN−q +

N ∑ q=0

)]

gq pbN+−1q

(

=−

1 2c

pbN+1 − pbN

Δt

+

pbN+−11 − pbN−1

Δt

)

.

(A.2)

14

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

Δt

Defining r = c 𝛽 , pulling out the q = 0 terms from the summation, and performing some algebraic steps simplifies Eq. Δx (A.2) to

( pbN+1

=

pbN

1−r r+1

)

(

+

1

)

r+1

(

(pbN−1

−

pbN+−11 )

−

r r+1

) (∑ N

gq pbN+−1q

q=1

+

N ∑

) gq pbN−q

.

(A.3)

q=1

Appendix B. Analysis using the integer-order partial reflection BC This supplementary study used the integer-order partial reflection BC given by Eq. (21) rather than the fractional BC to perform a similar analysis to that in § 3.1. Substituting Eq. (18) into Eq. (21) and evaluating at x0 = 0 yields j𝜔𝛾 (1 + R) + c(−jk + jkR) = 0.

(B.1)

Simplifying this equation and noting k = 𝜔∕c produces the formula

𝛾=

1−R . 1+R

(B.2)

A plot of the coefficient 𝛾 for the quadratic tapered ABH (reflection coefficient given by the black, solid curve in Fig. 4) is plotted in Fig B. 10. Note that the value of 𝛾 is complex.

Fig. B 10 Plot of the complex-valued coefficient 𝛾 for the integer-order partial reflection BC.

Having obtained the coefficient 𝛾 , we now develop the FD method for the integer-order partial reflection BC. It is exactly the same as the procedure given in §3.1.2, except the finite difference formula for the fractional Mur BC (Eq. (28)) is replaced by the FD Mur formula for the integer-order partial reflection BC, which is pbN+1 = pbN +

r − 𝛾 b+1 (p − pbN ) r + 𝛾 N −1

(B.3)

where the subindex N represents the last spatial node and r=c

Δt . (Δx)

(B.4)

The FD simulation for the partial reflection BC is performed across a wide range of frequencies to confirm its validity for the entire frequency range. The results are shown in Fig B. 11 by plotting the reflection coefficients from the TM method and the FD time domain simulation. Just as in Fig. 8, the reflection coefficient for the FD simulation was calculated from the ratio of amplitude of the FFTs at the corresponding wavenumber for the incident and reflected waves. The curves in Fig B. 11a nearly lie on top of each other. The percent difference between the curves is plotted in Fig B. 11b.

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

15

Fig. B 11 a) Plot of the reflection coefficients from the TM method and the FD simulation using the integer-order partial reflection BC. b) Plot of the percent difference error between the curves of the reflection coefficients using the integer-order partial reflection BC.

The percent error is always less than 0.05%. Comparing Fig. 8b to Fig. B. 11b, we observe that, while the error using the fractional BC was still low, the percent error in using the integer-order partial reflection BC is an order of magnitude less. Unfortunately, we were unable to definitively conclude why the fractional BC had a larger error than the integer-order partial reflection BC despite a thorough examination of the possibilities of error sources. We predict that this difference in error arises from a slight error in the calculation of a complex-valued fractional derivative, similar to what was found in Refs. [34,39]. One could argue against using a fractional BC due to the lower error associated with the integer-order partial reflection BC. However, recall that one of the main focuses of this paper was to show how fractional order mathematics could be used for modeling ABH terminations and to stimulate future thinking for the application of this mathematical tool to acoustic and vibration analysis. Our results using the fractional order operator certainly exhibited the ability of fractional order operators to serve as a fairly accurate means to model acoustic black hole terminations in ducts. References [1] M. Mironov, Propagation of a flexural wave in a plate whose thickness decreases smoothly to zero in a finite interval, Sov. Phys. Acoust. 34 (3) (1988) 318–319. https://www.scopus.com/record/display.uri?eid2-s2.0-3042686028origininward. [2] V.V. Krylov, Acoustic black holes: recent developments in the theory and applications, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 61 (8) (2014) 1296–1306, https://doi.org/10.1109/TUFFC.2014.3036, https://ieeexplore.ieee.org/document/6863850?arnumber6863850. [3] V. Krylov, F. Tilman, Acoustic black holes for flexural waves as effective vibration dampers, J. Sound Vib. 274 (3) (2004) 605–619, https://doi.org/10.1016/ j.jsv.2003.05.010, http://www.sciencedirect.com/science/article/pii/S0022460X03009490. [4] D. O’Boy, V. Krylov, V. Kralovic, Damping of flexural vibrations in rectangular plates using the acoustic black hole effect, J. Sound Vib. 329 (22) (2010) 4672–4688, https://doi.org/10.1016/j.jsv.2010.05.019, http://www.sciencedirect.com/science/article/pii/S0022460X10003457. [5] D. O’Boy, V. Krylov, Damping of flexural vibrations in circular plates with tapered central holes, J. Sound Vib. 330 (10) (2011) 2220–2236, https://doi.org/ 10.1016/j.jsv.2010.11.017. dynamics of Vibro-Impact Systems http://www.sciencedirect.com/science/article/pii/S0022460X10007790. [6] S. Conlon, J. Fahnline, F. Semperlotti, Numerical analysis of the vibroacoustic properties of plates with embedded grids of acoustic black holes, J. Acoust. Soc. Am. 137 (2015) 447–457, https://doi.org/10.1121/1.4904501. [7] H. Zhu, F. Semperlotti, Phononic thin plates with embedded acoustic black holes, Phys. Rev. B 91 (2015) 104304, https://doi.org/10.1103/PhysRevB.91. 104304, https://link.aps.org/doi/10.1103/PhysRevB.91.104304. [8] L. Zhao, F. Semperlotti, Embedded acoustic black holes for semi-passive broadband vibration attenuation in thin-walled structures, J. Sound Vib. 388 (2017) 42–52, https://doi.org/10.1016/j.jsv.2016.10.029, http://www.sciencedirect.com/science/article/pii/S0022460X16305648. [9] M.A. Mironov, V.V. Pislyakov, One-dimensional acoustic waves in retarding structures with propagation velocity tending to zero, Acoust Phys. 48 (3) (2002) 347–352, https://doi.org/10.1134/1.1478121. [10] O. Guasch, M. Arnela, P. Snchez-Martn, Transfer matrices to characterize linear and quadratic acoustic black holes in duct terminations, J. Sound Vib. 395 (2017) 65–79, https://doi.org/10.1016/j.jsv.2017.02.007, http://www.sciencedirect.com/science/article/pii/S0022460X17301013. [11] A.E. Ouahabi, D. O’Boy, V. Krylov, Experimental investigation of the acoustic black hole for sound absorption in air, in: Proceedings of the 22nd International Congress on Sound and Vibration, Florence, Italy, 2015. [12] A.E. Ouahabi, D. O’Boy, V. Krylov, Investigation of the acoustic black hole termination for sound waves in cylindrical waveguides, in: INTER-NOISE and NOISE-CON Congress and Conference Proceedings, Institute of Noise Control Engineering, San Francisco, USA, 2015. [13] K. Walther, Reflection factor of gradual-transition absorbers for electromagnetic and acoustic waves, IRE Transactions on Antennas and Propagation 8 (6) (1960) 608–621, https://doi.org/10.1109/TAP.1960.1144901. [14] V. Denis, A. Pelat, C. Touz, F. Gautier, Improvement of the acoustic black hole effect by using energy transfer due to geometric nonlinearity, Int. J. Non-Linear Mech. 94 (2017) 134–145, https://doi.org/10.1016/j.ijnonlinmec.2016.11.012. a Conspectus of Nonlinear Mechanics: A Tribute to the Oeuvres of Professors G. Rega and F. Vestroni http://www.sciencedirect.com/science/article/pii/S0020746216303675. [15] I. Podlubny, Fractional Differential Equations an Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications, Academic Press, New York, NY, 1999. [16] R. Herrmann, Fractional Calculus: an Introduction for Physicists, second ed., World Scientific Publishing Company, New Jersey, 2014, https://books.google. com/books?id60S7CgAAQBAJ. [17] K. Diethelm, The Analysis of Fractional Differential Equations, Springer, New York, NY, 2004.

16

J.P. Hollkamp and F. Semperlotti / Journal of Sound and Vibration 465 (2020) 115035

[18] V. Tarasova, V. Tarasov, Fractional dynamics of natural growth and memory effect in economics, European Research 23 (2016) 30–37, https://doi.org/10. 20861/2410-2873-2016-23-004. [19] H. Wang, J. Li, Surpassing the fractional derivative: concept of the memory-dependent derivative, Comput. Math. Appl. 62 (3) (2011) 1562–1567, https:// doi.org/10.1016/j.camwa.2011.04.028. [20] B. Narahari Achar, J. Hanneken, T. Clarke, Damping characteristics of a fractional oscillator, Physica A 339 (34) (2004) 311–319, https://doi.org/10.1016/j. physa.2004.03.030. [21] Y. Ryabov, A. Puzenko, Damped oscillations in view of the fractional oscillator equation, Physical Review B (Condensed Matter and Materials Physics) 66 (18) (2002), https://doi.org/10.1103/PhysRevB.66.184201 184201 1. [22] A. Tofighi, The intrinsic damping of the fractional oscillator, Physica A 329 (12) (2003) 29–34, https://doi.org/10.1016/S0378-4371(03)00598-3. [23] M. Meerschaert, R. McGough, Attenuated fractional wave equations with anisotropy, J. Vib. Acoust.. 136 (5). URL http://doi.org/10.1115/1.4025940. [24] P. Torvik, R. Bagley, On the appearance of the fractional derivative in the behavior of real materials, J. Applied Mechanics, Transactions ASME 51 (2) (1984) 294–298. [25] A.W. Wharmby, R. Bagley, Generalization of a theoretical basis for the application of fractional calculus to viscoelasticity, J. Rheol. 57 (5) (2013) 1429–1440. [26] B. Narahari Achar, J. Hanneken, Microscopic formulation of fractional calculus theory of viscoelasticity based on lattice dynamics, Phys. Scr. Volume T 2009 (T136) (2009) 014011–014018. [27] Z. Fellah, M. Fellah, C. Depollier, Transient wave propagation in inhomogeneous porous materials: application of fractional derivatives, Signal Process. 86 (10) (2006) 2658–2667, https://doi.org/10.1016/j.sigpro.2006.02.014. [28] G. Casasanta, R. Garra, Fractional calculus approach to the acoustic wave propagation with space-dependent sound speed, Signal, Image and Video Processing 6 (3) (2012) 389–392. [29] V. Tarasov, Acoustic waves in fractal media: non-integer dimensional spaces approach, Wave Motion 63 (2016) 18–22, https://doi.org/10.1016/j.wavemoti. 2016.01.003. [30] A. Carpinteri, F. Mainardi, Fractals and Fractional Calculus in Continuum Mechanics, Springer, London, 1997. [31] A. Carpinteri, B. Chiaia, P. Cornetti, The elastic problem for fractal media: basic theory and finite element formulation, Comput. Struct. 82 (6) (2004) 499–508, https://doi.org/10.1016/j.compstruc.2003.10.014. [32] J. Li, M. Ostoja-Starzewski, Fractal solids, product measures and fractional wave equations, in: Proceedings of the Royal Society: Mathematical, Physical and Engineering Sciences, vol. 465, 2009, pp. 2521–2536 (2108) http://rspa.royalsocietypublishing.org/content/465/2108/2521. [33] L. Cremer, M. Heckl, B. Petersson, Structure-Borne Sound: Structural Vibrations and Sound Radiation at Audio Frequencies, Springer, New York, 2005, https://books.google.com/books?idpAtwmxs-jaMC. [34] J.P. Hollkamp, M. Sen, F. Semperlotti, Analysis of dispersion and propagation properties in a periodic rod using a space-fractional wave equation, J. Sound Vib. 441 (2019) 204–220, https://doi.org/10.1016/j.jsv.2018.10.051, http://www.sciencedirect.com/science/article/pii/S0022460X18307314. [35] M. Israeli, S.A. Orszag, Approximation of radiation boundary conditions, J. Comput. Phys. 41 (1) (1981) 115–135. [36] J.R. Dea, Absorbing boundary conditions for the fractional wave equation, Appl. Math. Comput. 219 (18) (2013) 9810–9820, https://doi.org/10.1016/j.amc. 2013.03.113, http://www.sciencedirect.com/science/article/pii/S0096300313003743. [37] T.M. Atanackovi, S. Konjik, S. Pilipovi, D. Zorica, Complex order fractional derivatives in viscoelasticity, Mech. Time-Dependent Mater. 20 (2) (2016) 175–195, https://doi.org/10.1007/s11043-016-9290-3. [38] N. Makris, M. Constantinou, Models of viscoelasticity with complex-order derivatives, J. Eng. Mech. 119 (7) (1993) 1453–1464. [39] J.P. Hollkamp, M. Sen, F. Semperlotti, Model-order reduction of lumped parameter systems via fractional calculus, J. Sound Vib. 419 (2018) 526–543, https://doi.org/10.1016/j.jsv.2018.01.011, http://www.sciencedirect.com/science/article/pii/S0022460X18300117. [40] K. ichiro Hamanaka, Open, partial reflection and incident-absorbing boundary conditions in wave analysis with a boundary integral method, Coast Eng. 30 (3) (1997) 281–298, https://doi.org/10.1016/S0378-3839(96)00049-X, http://www.sciencedirect.com/science/article/pii/S037838399600049X. [41] G. Mur, Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic-field equations, IEEE Trans. Electromagn Compat. 23 (4) (1981) 377–382, https://doi.org/10.1109/TEMC.1981.303970. [42] E. Sousa, How to approximate the fractional derivative of order 1 < 𝛼 < 2, Int. J. Bifurc. Chaos Appl. Sci. Eng. 22 (4) (2012), https://doi.org/10.1142/ S0218127412500757 12500751250075.