Spectral analysis of acoustic plane waves in ductworks with uniform and varying cross-section waveguides

Article Journal of Vibration and Control 2024, Vol. 30(9-10) 2124–2138 © The Author(s) 2023 Article reuse guidelines: sagepub.com/journals-permissions DOI: 10.1177/10775463231175258 journals.sagepub.com/home/jvc Spectral analysis of acoustic plane waves in ductworks with uniform and varying cross-section waveguides Khaled M Ahmida Abstract The problem of low-frequency sound propagation in ducts involving varying cross-sections is analyzed. Spectral finite element formulations for the dynamic analysis of varying cross-section ducts are presented as tapered spectral elements. These elements can be used in connected waveguides through dynamic stiffness relations. Linear and polynomial duct types are considered. The elements’shape functions are established as well as the dynamic stiffness matrices describing the spectral relation between the acoustic pressure and the volume velocity. The presented spectral approach describes the acoustic wave motion at any position along arbitrary multiply connected duct waveguides within the one-dimensional planewave assumption. The derived expressions are verified experimentally, and by comparison with 3-D finite element solutions. The modeling of varying cross-section duct waveguides can play an important role in the design of acoustic excitation devices and horn-type loudspeakers. The advantage of the presented formulations lies in the low computational cost when compared to the finite element approach. Keywords Acoustic Ductwork, spectral elements, finite elements, wave propagation 1. Introduction One of the well-known ways to efficiently analyze vibrations of acoustic ducts with different boundary conditions and discontinuities is the use of matrix formulations such as the Finite Element Analysis (FEA). However, a large number of finite elements may be needed to adequately model a given dynamic system. Furthermore, the higher the frequency, the higher the number of elements needed for an acceptable solution. Analytical expressions can be found in the literature for the solution of the acoustic pressure in ducts with uniform cross-sections (Morse and Ingard, 1986;Kinsler et al., 1982;Nayfeh and Telionis, 1973). These uniform ducts can be easily modeled using spectral elements method (SEM); as well as finite elements. One class of problems of ducts is the case of the tubular and rectangular duct silencers used for the reduction of air-borne machinery noise and vent silencers. Hence, for the effective modeling of these types of ducts, it is necessary to formulate a tapered element that can easily be assembled with uniform cross-section elements. One known equation for solving the problem of low-frequency sound propagation in slowly varying ducts is Webster’s equation. This equation models the propagation of pressure waves in a hard-walled horn assuming that no transverse modes exist (Benade and Jansson, 1974). The main difference between the SEM and FEA is that in the spectral elements approach the calculations can only be carried out in the frequency domain, that is, there is no stiffness and mass matrices written separately, as in FEA, but rather a dynamic matrix. The tapered type of duct spectral elements presented in this paper deal with a symmetric variation of cross-sections. The advantage of formulating such elements is the capability of analyzing a multiply-connected network of acoustic ducts as connected elements in a relatively computationally cheap manner when compared to the FEA. Once the elementary dynamic matrix is obtained, the global matrix of the ducts network is built and solved for every frequency component. The disadvantage is that this kind of Department of Mechanical and Industrial Engineering, Faculty of Engineering, University of Tripoli, Tripoli, Libya Received: 20 November 2022; revised: 8 April 2023; accepted: 24 April 2023 Corresponding author: Khaled M Ahmida, Department of Mechanical and Industrial Engineering, Faculty of Engineering, University of Tripoli, University Campos, Tripoli, Libya. Email: [email protected] element describes the one-dimensional sound wave propagation while the FEA 3-D acoustic elements describe the acoustic pressure field in all directions inside a duct. This can be overcome by formulating a 2-D spectral element using analytical transverse modes as propagation modes as it is done, for example, with simply supported plate spectral elements [Lee and Lee, 1999;Doyle, 1997]. The limiting factor in the 2-D case is that there is a restricted number of boundary condition types that could be imposed on one of the sides of the element. On the other hand, the number of degrees of freedom used in FEA models could be reduced using a wave-based modal solution with an approach of mode matching or point collocation, using the so-called hybrid numerical method (Kirby, 2008). This paper is an amendment of a previous work conducted on tapered duct elements (Ahmida et al., 2003; Ahmida and Ferreira, 2004). In what follows, the element for uniform ducts is presented and Webster’s equation for tapered cross-sections is reviewed. The solutions for Webster’s equation for linearly-tapered and polynomially tapered elements are given. Based on these solutions, the tapered elements are formulated. The obtained expressions are validated experimentally by measuring the acoustic pressure in two connected cone-shaped ducts excited with a loudspeaker at one end. Then these are validated using finite element models via the general-purpose FEA program ANSYS ® . Different models are built to demonstrate the capability and robustness of the presented formulations. A numerical experiment is conducted to demonstrate the simplicity and validity of the given formulations through the dynamic analysis of a 3-D network of acoustic ducts, with a mix of uniform and tapered duct profiles. 2. Duct spectral element with uniform cross-section The uniform one-dimensional spectral duct element describes the acoustic wave propagation in ducts with a constant cross-sectional area. This formulation uses accurate shape functions based on Green’s function of the waveguide. The plane-wave propagation is described by Helmholtz and Euler equations, given by, respectively ∂2 ∂x2b Pðx,ωÞþk2b Pðx,ωÞ¼0 (1) ∂ ∂xPðx,tÞþρ∂u ∂t¼0 (2) with Pthe acoustic pressure, uthe particle velocity, ρthe mass density and k¼ω cthe wave number, cthe speed of sound and the symbol ⋀ representing frequency-dependent variables. The solution for the acoustic pressure and volume velocity is b Pðx,ωÞ¼A1eikx þB1eikx (3) b uðx,ωÞ¼ S ickρ ∂b PðxÞ ∂x(4) with Sthe cross-section area of the duct, i¼ffiffiffiffiffiffiffi 1 p, and A1& B1as constants obtained by applying the boundary conditions. For a spectral element of length L, consisting of two nodes and two DOF, the dynamic relation between the nodal acoustic pressure and the nodal volume velocity is then given by b u1 b u2¼S ρcðei2kL 1Þ1þei2kL 2eikL 2eikL 1þei2kL b P1 b P2(5) where b uðxÞ¼0atx¼0 and at x¼L, which represents a duct closed at both ends. Using the shape functions b g1UðxÞand b g2UðxÞ, the solution at any point xalong the duct is b Pðx,ωÞ¼b P1b g1UðxÞþb P2b g2UðxÞ(6) where the shape functions for a uniform duct element are given as exponential functions written simply as b g1Uðx,ωÞ¼eikx eikð2LxÞΔ b g2Uðx,ωÞ¼eikðLþxÞþeikðLxÞΔ(7) with Δ¼1ei2kL. These functions have a sinusoidal behavior, as shown in Figures 1 and 2for real and complex wave numbers, respectively. The wave number becomes a complex quantity when dissipation exists. In acoustics, the viscosity, the damping on the duct walls, and the fact that the process is adiabatic introduce dissipation terms in the equations describing wave propagation. The acoustic wave propagation process, based on experimental observations, is shown to be nearly adiabatic, as only insignificant amounts of energy can be exchanged between the fluid particles (Kinsler et al., 1982). The constant-amplitude behavior of b g1Uðx,ωÞand b g2Uðx,ωÞare observed in Figure 1. The imaginary part of the shape function has different behavior for complex wave numbers, as shown in Figure 2. 3. Tapered duct spectral element Low-frequency sound waves propagating in tapered ducts can be described by Webster’s equation. This is an ordinary one-dimensional differential equation that assumes plane wave approximation in hard-walled ducts, that is, the hypothesis of no transverse modes is adopted. This hypothesis is valid for the case where the corresponding wavelength is of the same order of magnitude as the transverse dimension of the duct. Considering the area variation S(x) along the Ahmida 2125 duct axis, Webster’s equation can be written as [Rienstra, 2002] ∂2P ∂t2¼c21 SðxÞ ∂ ∂xSðxÞ∂P ∂x (8) Using a time-harmonic standing waves solution for the pressure field, the Helmholtz equation for this tapered variation of S(x) is given by ∂ ∂xSðxÞ∂ ∂xb Pðx,ωÞþk2SðxÞb Pðx,ωÞ¼0 (9) In the following sections, two types of formulations are analyzed, the linear rectangular type and the polynomial type. 3.1. Cross-sections varying linearly Consider a linear variation of the cross-section area S of a duct of length L (Figure 3). The tapered duct has a constant depth, and it is connected to two other uniform ducts. The area change along the x-axis could be expressed by a linear function as SðxÞ¼m1ðxþε1Þ(10) Using the boundary conditions of S=S 0 at x= 0 and S= S L at x=L, the constants m1and ε1are found as ε1¼S0L SLS0 ,m1¼S0 ε1 (11) Figure 2. Shape functions for duct element with a uniform cross-section, and complex wave number (k=10-i0.1). Figure 1. Shape functions for duct element with a uniform cross-section, and real wave number k. 2126 Journal of Vibration and Control 30(9-10) The equation describing wave propagation for this kind of duct, using equations (9)and(10), is given by m1 d dxb PðxÞþm1ðxþε1Þd2 dx2b PðxÞþm1ðxþε1Þk2b PðxÞ¼0 (12) which has a Bessel-type solution given by b Pðx,ωÞ¼A2Y0ðkðxþε1ÞÞþB2J0ðkðxþε1ÞÞ (13) with J0the Bessel function of the first kind of zero-order and Y0the Bessel function of the second kind of zero-order. Using the conditions at the boundaries for the duct, that is, b P¼b P1at x=0andb P¼b P2at x = L, the constants A2and B2are determined. Rearranging the nodal variables, the dynamic stiffness matrix for a single duct element is then obtained as b u1 b u2¼hb KLib P1 b P2(14) where b uðxÞ¼0atx¼0 and at x¼L, which represents a duct closed at both ends. The explicit form of the [2 × 2] dynamic matrix ½b KLis given in Appendix A. The shape functions of this element are derived in the same manner as done before for the duct element with a uniform cross-section, and are given by These interpolation functions are real quantities when wave number kis real and become complex when damping is introduced, which represents the loss factor. The behavior of these two functions with respect to frequency is shown in Figure 4 for S 0 =10 4 m 2 , and S L =4×10 3 m 2 . When there are losses in the duct, the wave number becomes a complex quantity. Figure 5 shows the behavior of the two interpolation functions for k= 10-i0.1 m 1 . The effect of the tapered duct cross-section can be observed in the varying amplitude of the interpolation functions b g1Lðx,ωÞand b g2Lðx,ωÞ. For a real wave number, the real part has a decreasing pattern, clearly affected by the increasing cross-section area. This effect of damping is known as geometrical damping, where the decreasing amplitude is related to the fact that S 0 <S L . For a complex wave number, the real part maintains the same behavior but the complex part has a changing amplitude pattern. In both cases, the real part describes the acoustic wave propagation along the duct. 3.2. Cross-section variation of polynomial type This element is formulated in similar lines as done previously, with the area change described by the n-order polynomial SðxÞ¼mnðxþεnÞn(16) For the general case, that is, when n≥2, the corresponding variables and matrix formulations are given in Appendix B. Nonetheless, the special case of a second-order polynomial, that is, for n= 2, where the solution is simpler, is considered in what follows. Again, using the boundary conditions, the constants m2and ε2for the case of n= 2 are given by ε2¼1 22S0þ2ffiffiffiffiffiffiffiffiffi S0SL pL SLS0 ,m2¼S0 ε2 2 (17) Substituting into equation (9) yields 2m2ðxþε2Þd dxb PðxÞþm2ðxþε2Þ2d2 dx2b PðxÞ þm2ðxþε2Þ2k2b PðxÞ¼0 (18) which has a solution that could be written in the form b Pðx,ωÞ¼ 1 ðxþε2Þffiffiffi k pðA3cosðkðxþε2ÞÞþB3sinðkðxþε2ÞÞÞ (19) Figure 3. Acoustic duct with a linear variation of cross-section (shaded). b g1Lðx,ωÞ¼J0ðkðLþε1ÞÞ Y0ðkðxþε1ÞÞþY0ðkðLþε1ÞÞJ0ðkðxþε1ÞÞ Y0ðkðLþε1ÞÞJ0ðkε1ÞJ0ðkðLþε1ÞÞY0ðkε1Þ b g2Lðx,ωÞ¼ J0ðkε1ÞY0ðkðxþε1ÞÞY0ðkε1ÞJ0ðkðxþε1ÞÞ Y0ðkðLþε1ÞÞ J0ðkε1ÞJ0ðkðLþε1ÞÞ Y0ðkε1Þ (15) Ahmida 2127 Again, using the boundary conditions, A3and B3are determined and the dynamic relation is obtained as b u1 b u2¼hb KPib P1 b P2(20) where b uðxÞ¼0atx¼0andatx¼Lfor a duct closed at both ends. The explicit form of ½b KPis given in Appendix A. The shape functions are derived as before, and are given by Figure 5. Shape functions for a tapered element with a linear variation of the cross-section area, and complex wave number k. b g1Pðx,ωÞ¼ε2ðcosðkðxþε2ÞÞsinðkLÞ cosðkε2Þþcosðkðxþε2ÞÞ cosðkLÞ sinðkε2ÞÞ ðxþε2Þ sinðkLÞ ε2ðsinðkðxþε2ÞÞ cosðkLÞ cosðkε2Þþsinðkðxþε2ÞÞ sinðkLÞ sinðkε2ÞÞ ðxþε2Þ sinðkLÞ b g2Pðx,ωÞ¼ðLþε2Þðsinðkε2Þ cosðkðxþε2ÞÞþcosðkε2Þ sinðkðxþε2ÞÞÞ ðxþε2Þ sinðkLÞ (21) Figure 4. Shape functions for a tapered element with a linear variation of the cross-section area, and real wave number k. 2128 Journal of Vibration and Control 30(9-10) The interpolation functions b g1Pðx,ωÞand b g2Pðx,ωÞare real functions for real wave numbers and complex when kis complex. The characteristics of these functions are similar in behavior to those with linear variations of cross-sections (Figures 4 and 5), demonstrating the effect of geometrical damping. Large numbers of cross-section variations can be represented using equation (16) with higher-order polynomials, that is, with n> 2. The solution would then be of Bessel-type for odd values of nand a complex exponential type for even values of n. Other types of function representations can also be considered, for instance, the exponential function. Nevertheless, the used formulations cover a wide range of flare configurations. For steep cross-section areas, transverse modes could appear across the duct axis at lower frequencies. These modes are not modeled by the spectral elements developed in this paper and new formulations need to be developed to account for the transverse waves in ducts. 4. Experimental validation This experiment aims to validate the obtained expressions for a linear type of tapered duct. Two plastic conic-shaped acoustic ducts were connected as shown in Figure 6. These are similar having a wall thickness of 2 mm and cross-section areas of Sð1Þ 0=0.0227m 2 ,Sð2Þ 0¼Sð1Þ L=0.0028m 2 and Sð2Þ L=0.0487m 2 . A loudspeaker was used as an excitation, fitattheendofareaSð1Þ 0. Two small holes were made on the system surface for the two ICP ® microphones (PCB model TMS130A10) to measure the acoustic pressure. Point #1 is located at the junction point and the other is located in the middle of cone #2 (Figure 6(a)). The cut-off frequency, above which the sound waves propagate not only in the longitudinal direction but also in the transverse direction, was calculated. In a circular-crosssection duct, the cut-off frequency in Hz is fcutoff ¼0:586 c=D, where c is the speed of sound in m/s, and Dis the largest diameter of the waveguide in m (Kinsler et al., 1982). In our case, the cut-off frequency was found equal to 800 Hz. Accordingly, a frequency range of 1 Hz to 800 Hz was used for the analysis. An HP data acquisition system was used to acquire the measured signals. The pressure frequency response was measured at points #1 and #2. The reference signal was the velocity of the loudspeaker cone measured with a Laser Doppler Vibrometer (Figure 6(b)). The two-duct system was also modeled via SEM using the aforementioned expressions for the linearly tapered element. The numerical model consists of only two Figure 6. The experimental setup of two connected tapered ducts: (a) experimental setup; (b) schematic diagram. Ahmida 2129 elements with only three DOFs to solve. It should be mentioned that the open-ended condition is not the zero impedance condition. At this end, the acoustic radiation impedance Zoshould be added since the duct is open and radiates sound into the surrounding medium. For an unflanged open duct, experiments and theory indicates that the radiation impedance could be approximated by (Kinsler et al., 1982) Zo¼ρcS1 4k2R2þi0:6kR(22) where Ris the radius of the open end, kis the wave number, cis the speed of sound, Sis the cross-section area of the open end, and ρis the air density. This equation is consistent with the assumption that the wavelength is large compared to the largest radius of the duct system. The frequency responses were calculated at points #1, and at point #2 using the shape functions. These are shown in comparison with the measured ones in Figures 7 and 8. The good agreement observed between the measured and numerically modeled frequency response at point #1 was used for the verification of the stiffness matrices for the tapered element and the agreement observed at point #2 was used to verify the obtained shape functions. The small difference observed in the amplitudes can be related to the constant value of the loss factor used in the SEM model. In this case, a loss factor of 0:03 was used for the calculation of wave number k. 5. Validation via finite element analysis The aforementioned dynamic relations are verified by using 3-D finite element solutions and are judged by comparison of the natural frequencies of the duct element. The two types of tapered ducts were modeled using the commercial FEA program ANSYS ® . Slow and steep varying cross-sections were considered. The FEA model assumes the Neumann boundary condition of hard-walled ducts closed at both ends, as is the case for the spectral element model. It is important to mention that in all the results presented for linear and polynomial tapers, only one single spectral element was used to represent the tapered duct. A finite element acoustic modal analysis was conducted to calculate the first few natural frequencies and mode shapes. All the existing modes in the analyzed frequency range are found; whether longitudinal or transverse modes. In all following analyses, the medium is air having a mass density of 1.21 kg/m 3 , and a speed of sound of 340 m/s. 5.1. Linear-type waveguide model Consider a tapered waveguide of length L= 2m with a cross-section area varying linearly and symmetrically along its axis. Two kinds of duct shapes are analyzed: slowly tapered and steeply tapered. The FEA element used is the FLUID220 higher order 3-D 20-node solid element, used for modeling the fluid medium. 5.1.1. Slow taper. The cross-section areas at the ends of the duct have a constant depth of 0.1 m and a width varying from 0.1 m at one end to 0.3 m at the other end, thus resulting in the cross-section areas of S 0 = 0.01 m 2 at x= 0 and S L = 0.03 m 2 at x=L. A modal analysis was conducted and it was found that the first few natural frequencies converged to an FEA model of 26,493 elements with an element size of 0.02 m. The identified natural frequencies are listed in Table 1, compared with the SEM ones. The first natural frequency is zero as the first mode shape describes the constant pressure mode, which is equivalent to a rigid body mode in structures. Figure 7. Frequency response of acoustic pressure at point #1: via SEM (straight line); measured (dashed line). Figure 8. Frequency response of acoustic pressure at point #2: via SEM (straight line); measured (dashed line). 2130 Journal of Vibration and Control 30(9-10) In the SEM model, a constant-frequency spectrum was used as an excitation applied at the section of the smaller area of the duct. This excitation is a harmonic velocity (piston-like) with unit magnitude. Observe that the natural frequency values of the spectral element and FEA/ANSYS ® compare very well. The SEM frequency response function was calculated through the direct inverse method. This is numerically cheap as the matrices involved have dimensions of only 2 × 2. Then the natural frequencies were identified using a specific code written in MATLAB ® , see Figure 9. It is important to mention that the eighth mode shape, according to FEA/ANSYS model is at 611.3 Hz, and it is the first non-plane mode. This mode could not be identified by the presented SEM formulation as it has a non-plane wave propagation pattern. This mode shape is characterized by the transversal propagation normal to the duct section, in contrast to the third mode shape which is an axial mode, as shown in Figure 10. 5.1.2. Steep taper. The steep taper used in this analysis has a constant depth of 0.2 m and the width varies from 0.2 m at one end to 1.0 m at the other end, hence, resulting in the cross-section areas of S 0 = 0.04 m 2 at x= 0 and S L = 0.2 m 2 at x=L. The first few natural frequencies converged to an FEA model of 18,180 elements. The acoustic modal analysis was performed and the identified natural frequencies are listed in Table 2 for comparison with the SEM model. Only one spectral element is needed to model this tapered duct and obtain the natural frequencies. Due to the noticeable difference in the cross-sectional areas of the duct at the two ends, wherein this case S L is 5 times S 0 , the FEA model calculated some transverse modes. The spectral element model does not predict these transverse modes shown in Figure 11(b). These transverse modes emerged because, in this case, at a certain distance along the taper the width of the duct has a considerable increase, with respect to its height. This increase affects the cut-off frequency above which these transverse modes start to appear. In this case, it is only ∼200 Hz. It is important to highlight that as long as the plane wave approximation is valid, the validity of SEM solutions is not affected and accurate results can still be obtained, because it is based on exact solutions of the wave equation, see error percentages in Table 2. Nevertheless, one issue to question and needs extensive investigation is how slow should the Table 1. Natural frequencies identified for the linear-type slowly tapered duct. Mode no. via SEM, Hz via FEA/ANSYS®, Hz 1 88.51 88.47 2 172.02 171.94 3 256.39 256.28 4 341.06 340.91 5 425.85 425.67 Figure 9. The FRF calculated via SEM at the excitation node of a linear-type tapered duct. Ahmida 2131 variation of the duct cross-section be for valid approximations. In our case of the slow taper above, where the ratio of cross-sectional areas is 3, the first non-plane mode shape observed was the eighth mode shape (611.3 Hz, not shown in Table 1). Our formulations are based on Webster’s equation, hence the issue of how slow the taper should be perhaps requires further investigation of the limitations of Webster’s equation itself, with respect to abrupt changes in the cross-section, or when the acoustic disturbance takes the form of a non-plane wave. Experimentally supported research would probably be needed for that matter. 5.2. Polynomially varying cross-section waveguide model In this case, the duct cross-section area varies according to a second-order polynomial. Only one steep tapered duct waveguide of length L= 1m is analyzed. The duct element has a circular cross-section area that varies symmetrically along its axis. The steep taper used has a cross-section area of S 0 = 0.001 m 2 at x= 0 and a larger area of S L = 0.01 m 2 at x=L. The first few modes converged using an FEA model of 8112 elements with an element size of 0.02 m. These frequencies are given in Table 3, compared to the oneelement SEM model. Observe that the natural frequencies are well-identified using just one spectral element, instead of solving a relatively large FEA matrix system. The differences in the natural frequencies obtained are very low for the first few modes; about 2%. This tiny discrepancy could be related to the SEM solution that is based on Bessel functions, and to Figure 10. Two distinct mode shapes of the FEA model: (a) the third mode shape (plane mode); (b) the eighth mode shape (non-plane mode). The duct is linear-type slowly tapered. DOF is the acoustic pressure. Table 2. Natural frequencies identified for the linear-type steeply tapered duct. Mode no. via SEM, Hz via FEA/ANSYS®, Hz Difference 1 91.68 91.1 0.63% 2 174.36 173.1 0.72% 3 Non-plane 201.9 — 4 258.15 256.3 0.71% 5 Non-plane 300.1 — 6 342.45 339.5 0.86% 7 Non-plane 378.7 — 8 Non-plane 383.1 — 9 426.99 424.6 0.56% 2132 Journal of Vibration and Control 30(9-10)