- Email: [email protected]
Boundary element analysis of bar silencers using the scattering matrix with two-dimensional finite element modes
Engineering Analysis with Boundary Elements 74 (2017) 100–106
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
crossmark
Boundary element analysis of bar silencers using the scattering matrix with two-dimensional finite element modes ⁎
L. Yanga, P. Wangb, T.W. Wub, a b
School of Power and Energy Engineering, Harbin Engineering University, Harbin, Heilongjiang 150001, PR China Department of Mechanical Engineering, University of Kentucky, Lexington, KY 40506-0503, USA
A R T I C L E I N F O
A BS T RAC T
Keywords: Acoustics Mufflers Silencers Boundary element method Finite element method
Bar silencers used in industry may consist of a large array of rectangular or round bars packed in a rectangular lattice arrangement. Due to the size of the lattice, normally only a single unit that represents a building block for the lattice is isolated for analysis purposes. Even with one isolated unit, the inlet and the outlet are still quite large, and the plane-wave cutoff frequency can be very low. Therefore, higher-order modes must be considered at the inlet and outlet in order to calculate the transmission loss. This paper uses the recently developed “impedance-to-scattering matrix method” to convert the element-based impedance matrix into the mode-based scattering matrix for transmission loss calculation. Depending on the shape of the inlet and outlet, it may not always be possible to find an analytical expression of the modes needed for the modal expansion. In this paper, the two-dimensional finite element method is used to extract the eigenvalues and the eigenvectors of the inlet/ outlet cross section. The eigenvectors are then used in the modal expansion to convert the impedance matrix into the scattering matrix. Test cases include several commonly used inlet and outlet configurations, such as rectangular, circular and triangular cross sections.
1. Introduction Dissipative silencers are widely used to attenuate broadband noise generated from gas turbine engines, power plants and HVAC ducts. There are a variety of designs of dissipative silencers distributing sound absorbent in various ways. A very simple design called the lined duct design places the sound absorbing material around the circumference of the cross-sectional area. Another prevailing design is the splitter silencer or the parallel-baffle silencer, which arranges multiple sound absorbing baffles parallel to the direction of exhaust. Normally, the sound attenuation of such silencers is proportional to the perimeterarea ratio and length [1]. There have been plenty of research activities on lined-duct silencers or splitter silencers either analytically, numerically or experimentally [2–7]. In 1983, Nilsson and Söderqvist proposed the idea of bar silencers and claimed that an array of square bars made of sound absorbing materials have certain advantages over a similarly configured splitter silencer [8]. In 1996, Cummings and Astley [9] investigated the acoustical behavior of square bar silencers using the finite element method (FEM) and compared to the experimental data. To date, the FEM is still the most developed and widely used numerical technique for large dissipative silencers. Kirby and his co-
⁎
workers [10–14] used a hybrid analytical-FEM to study the acoustical performance of various large dissipative silencers. To apply the hybrid technique, the 2D FEM is first employed to extract the eigenvalues and the associated eigenvectors of an axially uniform cross section. These 2D transversal modes are then used in the modal expansion along the axial direction if the cross section remains the same. To determine the unknown amplitudes in the modal expansion, either a point collocation method or an integral-based mode matching scheme is adopted to enforce the continuity of sound pressure and particle velocity at both ends where the uniform section meets the flanges or any irregular junctions. Since the FEM is mainly used on a 2D cross section to extract the modes, the hybrid analytical-FEM is a very efficient numerical technique for silencers with a very long axially uniform section. As a viable alternative to the FEM, the boundary element method (BEM) has also been used for muffler and silencer analysis [15–17]. The direct mixed-body BEM in Refs. [15–17] can handle complex internal components, such as thin bodies, perforated tubes, and multiple bulk-reacting materials, in one single BEM domain without resorting to the multi-domain approach. It is well known that the BEM has the advantage of modeling the surface only. However, for a large silencer, even the surface mesh may still contain too many elements for
Corresponding author. E-mail address: [email protected] (T.W. Wu).
http://dx.doi.org/10.1016/j.enganabound.2016.11.001 Received 30 January 2016; Received in revised form 29 September 2016; Accepted 2 November 2016 0955-7997/ © 2016 Published by Elsevier Ltd.
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
a desktop computer to run comfortably. Normally, one way to solve large-scale problems in BEM is to use the fast multipole method (FMM) [18,19]. The FMM is especially useful for exterior radiation and scattering problems. However, the FMM has not been extended to silencer problems with complex internal components and multiple bulk-reacting materials yet. Fortunately, the interior domain of a large silencer can always be divided into several smaller substructures by using the substructure BEM [20] so that each substructure can fit within the memory limitation of a desktop computer. Each substructure produces an impedance matrix, and a so-called “impedance matrix synthesis” procedure [20] is used to combine all the substructures. Until very recently, the BEM analysis of mufflers and silencers has been limited to problems with a small inlet and a small outlet only due to its dependency on either the four-pole matrix or the anechoic termination boundary condition for transmission loss (TL) calculation [15,16]. Without a proper introduction of modal expansion into the BEM, it is difficult to model large silencers with many higher-order modes emerging at the inlet and outlet cross sections. Zhou et al. [21] recently proposed a reciprocal identity method in conjunction with the BEM impedance matrix to extract the higher-order modes at the inlet and outlet. Each reciprocal identity couples the analytical modal expansion in the inlet and outlet ducts to a BEM solution with a random boundary condition set. The BEM impedance matrix can naturally provide more than enough such solutions for the reciprocal identity coupling. The method can be regarded as an indirect postprocessing filter applied to the BEM impedance matrix to extract the higher-order modes for TL calculation. Parallel to the reciprocal identity method, a different technique called the “impedance-to-scattering matrix method” was recently developed by Wang and Wu [22] as a direct collocation approach to convert the element-based impedance matrix into the mode-based scattering matrix for TL computation. The BEM impedance matrix relates sound pressures at the inlet and outlet to the corresponding particle velocities, while the scattering matrix relates the modes at the inlet and outlet [23]. It should be noted that the term “scattering matrix” can be a little confusing and it is not the same as the T-matrix method for exterior scattering problems. The scattering matrix is a special tool designed for muffler and silencer analysis, and is equivalent to the four-pole matrix if the plane-wave assumption is still valid at the inlet and outlet. In the impedance-to-scattering matrix method [22], each sound pressure and particle velocity can be directly expanded in terms of the duct modes at the centroid of each constant boundary element. These point-wise expansions are then related by the BEM impedance matrix, and the scattering matrix can be obtained after a few matrix operations. Like the reciprocal identity method, the impedance-to-scattering matrix method can be regarded as a postprocessing filter applied to the BEM impedance matrix. The main difference between these two methods is that the reciprocal identity method is an integral-based method, while the impedance-to-scattering matrix method is a direct collocation method. Both methods rely on an analytical modal expansion at the inlet and outlet. Normally this is not an issue for a rectangular or circular inlet/outlet. However, more complicated inlet/outlet configurations will require a numerical solution of the duct modes first. In this paper, the impedance-to scattering matrix method by Wang and Wu [22] is expanded to include the use of numerical modes from the 2D FEM so that any irregular inlet/outlet configurations can be modeled. The function of the 2D FEM is to extract the eigenvalues and the associated eigenvectors (modes) of the inlet/outlet cross sections. The 2D FEM procedure is the same as the one used in the hybrid analytical-FEM [10–14] or the hybrid analytical-BEM [24], except that it is only done at the inlet/outlet cross sections as a post-processing tool to extract the duct modes. The proposed method is by no means a hybrid technique because the silencer structure itself is still modeled by the same substructure BEM as in Ref. [20] with the impedance matrix being the primary output from the BEM. Modal expansion and the
Fig. 1. Flowchart of the computation procedure.
associated 2D FEM analysis are only introduced in the end at the inlet and outlet cross sections after the resultant BEM impedance matrix is produced. Fig. 1 shows the flowchart of the proposed computation procedure. Figs. 2–4 show three typical lattice arrangements of bar silencers. Usually a small building block is isolated from the lattice for analysis purposes. The building block in Figs. 2 and 3 is simply a round or rectangular bar housed in a rectangular duct. The duct walls are assumed rigid due to symmetry. The building block in Fig. 4 actually begins with a hexagon duct. Due to rotational symmetry, the hexagon eventually reduces to just a triangle. Although rectangular and circular
Fig. 2. A rectangular module isolated from an aligned lattice arrangement of round bars.
101
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
subdomain integral equations are summed together, while Eq. (2) is the normal-derivative integral equation derived from Eq. (1).
⎞ ∂G 0 (rp, rQ ) + + jk 0 z 0 vn (rQ ) G0 (rp,rQ ) ⎟ dS + ∫ (p − p− ) dS ∂nQ T ⎠ ∼ ∼ ⎞ ⎛ ∂G (rp, rQ ) ∂G (rp, rQ ) + ∼ ∼ + ∫ ⎜p (rQ ) ∂n + jk z∼vn (rQ ) G (rp,rQ ) ⎟ dS + ∫ (p − p− ) dS ∂nQ B⎝ BTB Q ⎠ ∼ ∼ ⎡ ⎛ ∂G 0 (rp, rQ) ∂G (rp, rQ ) ⎞ ∂G (rp, rQ ) ∼ + ∫ ⎢p0 ⎜ ∂n − ∂n + jk 0 z 0 vn (rQ ) G0 (rp, rQ ) ⎟ − z 0 ζ vn (rQ ) ∂n IP ⎣ ⎝ Q Q Q ⎠ ⎛
∫R ⎜p (rQ ) ⎝
∂G 0 (rp, rQ ) ∂nQ
⎛ ∂G 0 (rp, rQ) ∼ ∂G∼(rp, rQ) ⎞ ∼ ∼ − jk z∼vn (rQ ) G (rp, rQ )] dS + ∫ ⎜p0 ∂n − p ∂n ⎟ dS ATB ⎝ Q Q ⎠ ⎧ p (rp ), P ∈ Ω. ⎪ ⎪ 0.5p (rp ), P ∈ R + B. ⎪ + − = ⎨ 0.5[p (rp ) + p (rp )], P ∈ T + BTB. ⎪ ∼ ⎪ p0 (rp ) + 0.5z 0 k 0 ζ vn (rp ), P ∈ IP. ⎪ 0.5[p (rp ) + ∼ p (rp )], P ∈ ATB. ⎩ 0
Fig. 3. A rectangular module isolated from an aligned lattice arrangement of square bars.
⎛
∂ 2G 0 (rp, rQ )
⎝
∂nQ ∂np
∫R ⎜p (rQ )
∂G 0 (rp,rQ ) ⎞
⎟ dS + ∫ T ⎠ ∼ ⎞ ∂G (rp,rQ ) ∼ + jk z∼vn (rQ ) ∂n ⎟ dS p ⎠
+ jk 0 z 0 vn (rQ )
∼ ⎛ ∂ 2G (rp, rQ ) + ∫ ⎜p (rQ ) ∂n ∂n Q p B⎝ ⎡ ⎛ ∂ 2G 0 (rp, rQ) + ∫ ⎢p0 ⎜ ∂n ∂n − Q p IP ⎣ ⎝
∂ 2G 0 (rp, rQ )
∂np
∂nQ ∂np
(1)
(p+ − p− ) dS
∼ ∂ 2G (rp, rQ ) ⎞ ∂nQ ∂np
⎟ ⎠
∼ ∼ ∂ 2G (rp, rQ ) ∂G 0 (rp, rQ ) ∂G (rp, rQ ) ⎤ ∼ ∼ − z 0 ζ vn (rQ ) ∂n ∂n + jk 0 z 0 vn (rQ ) ∂n −jk z∼vn (rQ ) ∂n ⎥ dS Q p p p ⎦
⎛ ∂ 2G 0 (rp, rQ) ∼ ∂ 2G∼(rp, rQ) ⎞ + ∫ ⎜p0 ∂n ∂n − p ∂n ∂n ⎟ dS + ∫ BTB ATB ⎝ Q p Q p ⎠ ⎧ 0, P ∈ T + ATB + BTB. =⎨ ∼∼ ⎩−0.5(jk 0 z 0 + jk z ) vn (rp ), P ∈ IP.
∼ ∂ 2G (rp, rQ ) ∂nQ ∂np
(p+ − p− ) dS
⎪ ⎪
(2) ∼ ∼ − jkr PQ where G0 = and G = e /4πrPQ are the free-space Green's functions in the air and the sound-absorbing material, respectively, P and Q are the collocation point and the integration point, respectively, nP and nQ are the corresponding normal vectors, rPQ = rP − rQ , z 0 and z∼ represent the characteristic impedances of air ∼ and the sound-absorbing material, k 0 and k the wavenumbers of ∼ corresponding media, respectively, ζ the transfer impedance of any perforated interfaces. The symbols R, T, B, IP, ATB and BTB represent “regular”, “thin”, “bulk-reacting”, “interface with a perforated facing sheet”, “air-thin-bulk reacting” and “bulk reacting-thin-bulk reacting” surfaces as defined in Fig. 5. p denotes sound pressure (on IP or ATB, ∼ are the sound pressure variables in the air and the soundp0 and p absorbing material, respectively) and vn is the normal particle velocity. The unit normal vector on either R or B is pointing into the interior acoustic domain. On T and BTB, the unit normal vector can be pointing into either side of the thin surface. The side into which the normal is pointing is called the positive side while the opposite side is called the negative side. The corresponding sound pressure values are denoted by p+ and p−, respectively. The unit normal vector on IP or ATB is pointing
e−jk 0 rPQ /4πrPQ
Fig. 4. A triangular module isolated from a shifted lattice arrangement of round bars.
duct modes are well known analytically, triangular or hexagon duct modes may not be easy to find. In some single-silencer applications, the inlet and/or outlet may also be an ellipse, trapezoid, or other irregular shapes. In this paper, the 2D FEM is used to extract the duct modes at the inlet and outlet, regardless of the shape. The method is first verified by comparing to an analytical solution for a single round bar silencer housed in a circular duct. The result of a rectangular bar-silencer test case is validated by comparing to an existing hybrid analytical-FEM solution available in literature. Then a triangular test case follows and is compared to the rectangular design.
2. Theory 2.1. Direct mixed-body BEM The concept of the direct mixed-body BEM [15–17,20] actually begins with the conventional multi-domain BEM by subdividing the acoustic domain into several well-defined homogeneous subdomains. The Helmholtz integral equation may be written for each individual subdomain. All the integral equations are then summed to create a single integral equation. The normal-derivative hypersingular integral equation is used to provide an additional equation at any interfaces that have two unknown variables. A review of hypersingular integrals can be found in Chen and Hong [25]. The direct mixed-body BEM is ideally suited to mufflers and silencers with complex internal components. The integral equations needed in this paper are listed in Eq. (1) and Eq. (2) and the complete edition can be found in Ref. [20]. Eq. (1) is the conventional Helmholtz integral equation except that all the
Fig. 5. Types of surfaces defined in the direct mixed-body BEM.
102
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
circular or rectangular inlet/outlet, the eigenfunctions can be obtained analytically. For irregular inlet/outlet shapes, the eigenvalues and the corresponding eigenfunctions must be obtained numerically. In the paper, the 2D FEM is used to calculate the eigenvalues and the associated eigenfunctions of the inlet/outlet cross sections. The governing differential equation for the transversal sound pressure in a 2D cross section is [26] 2 ∇2xy pxy + kxy pxy = 0
(6)
where pxy is the transversal sound pressure, and kxy is the transversal wavenumber that must satisfy the following equation: 2 kxy + kz2 = k 02
In FEM, the transversal sound pressure at any point is expressed in terms of the shape functions
Fig. 6. Substructure technique.
pxy = N Tp
into the air side.
(8)
where the matrix N consists of the column vectors of the global shape functions, and p is the column vector of the nodal values of transversal sound pressure. Applying the rigid-wall boundary condition and the Galerkin procedure [26] to Eq. (6) yields the following eigen-equation:
2.2. Impedance matrix and substructure BEM For large silencers, the substructure BEM [20] is a very useful technique as it can be difficult to analyze a large silencer in one single BEM model on a desktop computer. A demonstrated in Fig. 6, a large silencer can be divided into several smaller substructures so that each substructure can fit within the memory limitation of a desktop computer. Specifically, if a silencer has an axially uniform section, the impedance matrix of a small template taken from the uniform section can be repeatedly used without having to model the long uniform section. The impedance matrix of each substructure relates sound pressures at the inlet and outlet of the substructure to the corresponding particle velocities. Multiple computers may be used simultaneously for different substructures as each substructure analysis is independent of the others. All the substructure impedance matrices are then combined to form the resultant impedance matrix for the whole silencer. The detailed procedure is described in Ref. [20]. The end result is the resultant impedance matrix for the large silencer:
⎡ p I ⎤ ⎡ Z11 Z12 ⎤ ⎡ vI ⎤ ⎥⎢ ⎥ ⎢p ⎥ = ⎢ ⎣ O ⎦ ⎣ Z21 Z22 ⎦ ⎣ vO ⎦
(7)
2 (K − kxy M) p = 0,
(9)
(∇N)e (∇N)Te dSe Se
where K = ∑ ∫
(N)e (N)Te dSe Se
and M = ∑ ∫
are the
stiffness matrix and the mass matrix, respectively, subscript "e" denotes the element. If the number of nodes is n, n eigenvalues (kxy )i (1 ≤ i ≤ n ) and the associated eigenvectors (Φ)i with length of n may be obtained by solving Eq. (9). Assemble all the eigenvalues and eigenvectors in the vector k and the matrix Φ , respectively. Let A denote the wave amplitudes at the inlet and C the wave amplitudes at the outlet. Apply Eqs. (4) and (5) to the inlet and outlet individually, and sound pressures and particle velocities at the inlet and the outlet may be written in terms of their corresponding modes as
(3)
where the vectors p and v represent the sound pressures and particle velocities, and the subscripts I and O denote the inlet and the outlet, respectively.
ΦAn ⎤ ⎡ A+ ⎤ ⎡ p I ⎤ ⎡ ΦAn n ⎢ v ⎥ = ⎢ kAzn Φ − kAzn Φ ⎥ ⎢ − ⎥ ⎣ I ⎦ ⎢⎣ ρ ω An ⎥ ρ0 ω An ⎦ ⎣ An ⎦ 0
(10)
ΦCn ⎤ ⎡ C+ ⎤ ⎡ pO ⎤ ⎡ ΦCn ⎥⎢ n ⎥ k ⎢ v ⎥ = ⎢ kCzn Φ ΦCn ⎥ ⎣ C−n ⎦ ⎣ O ⎦ ⎢⎣ ρ ω Cn − ρCzn ⎦ 0 0ω
(11)
where kAzn and k Czn contain all the eigenvalues, and ΦAn and ΦCn contain the corresponding eigenvectors. From Eqs. (3), (10), and (11), the following relationship is obtained
2.3. Impedance-to-scattering matrix method with FEM modes
⎤ ⎧ C+n ⎫ ⎡ MPI−Z11MVI ⎤ ⎧ A+n ⎫ ⎡ Z12 MVO ⎥⎨ ⎬ ⎢ ⎥ ⎨ −⎬ = ⎢ ⎣ Z21MVI ⎦ ⎩ An ⎭ ⎣ MPO−Z22 MVO ⎦ ⎩ C−n ⎭
The impedance matrix in Eq. (3) can be converted into the scattering matrix by expressing each sound pressure and particle velocity at the inlet and outlet in terms of the duct modes [22]. In this paper, the duct modes are obtained by the 2D FEM so that the impedance-to-scattering matrix method is not limited to certain regular shapes only. At the inlet and outlet, sound pressure and particle velocity at any point nC can be expressed in a modal expansion:
(12)
where
MPI = [ ΦAn ΦAn ] MPO = [ ΦCn ΦCn ]
⎡k ΦAn − MVI = ⎢ ρAzn ⎣ 0ω
⎤ kAzn Φ ρ0 ω An ⎥ ⎦
M
pnC =
∑ ΦnC m (Fm+e−jkzm z + Fm−e jkzm z ) m =1
⎡k ΦCn − MVO = ⎢ ρCzn ⎣ 0ω
(4)
and
vnC =
1 ρ0 ω Fm+
Although Eq. (12) may look like a “transfer scattering matrix” that relates the wave amplitudes at the inlet to the wave amplitudes at the outlet, matrix inverse should not be performed at this stage because Eq. (12) does not represent a well-posed boundary value problem. A wellposed boundary value problem should have a known condition at the inlet and another known condition at the outlet. For example, to find the TL, a certain incident wave should be given at the inlet and the outlet is assumed anechoic. Rearrange Eq. (12) to get
M
∑ ΦnC m kzm (Fm+e−jkzm z − Fm−e jkzm z ) m =1
⎤ k Czn ΦCn ⎥ ρ0 ω ⎦
(5)
Fm−
where and are the modal amplitudes corresponding to the acoustic waves of order m travelling in the positive and negative z directions, respectively, kz represents the wavenumber in the axial direction, ΦnC m denotes the eigenfunction value at the point nC . For a 103
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
⎧ C+ ⎫ ⎧ A+ ⎫ [SO] ⎨ −n ⎬ = [SI] ⎨ −n ⎬ ⎩ An ⎭ ⎩ Cn ⎭
(13)
Since there are always more boundary elements at the inlet and outlet than the number of modes, a least-square matrix inverse may be performed on Eq. (13) to obtain the scattering matrix S :
⎧ C+n ⎫ ⎧ A+ ⎫ ⎨ − ⎬ = [S] ⎨ −n ⎬ ⎩ An ⎭ ⎩ Cn ⎭
(14)
2.4. Transmission loss To calculate the transmission loss, a single incident plane wave is prescribed at the inlet, and the outlet is assumed anechoic. In other words, let A+n = [1, 0, 0, ... ,0]T and C−n = [0, 0, 0, ... ,0]T . The transmitted wave amplitudes C+n and the reflected wave amplitudes A−n can be obtained from the scattering matrix relationship, Eq. (14). The transmission loss is defined as the difference of sound power in dB between the single incident plane wave and the transmitted waves. That is,
TL = 10 log10 (WI / WO )
Fig. 8. TL comparison of the round bar silencer.
It is noted that a similar test case was also used in Refs. [21,22] to validate the reciprocal identity method and the impedance-to-scattering matrix method, respectively. The purpose of repeating this particular test case here is to validate the proposed scattering matrix conversion method using the 2D FEM modes. The round bar is made of polyester with flow resistivity R=16000 Rayls/m, and is covered by a 30% open perforated facing sheet to protect the material from being blown away by the exhaust gas. Both ends of the bar are covered by a steel plate. With reference to Fig. 6, the dimensions used in this test case are r1 = 0.342m , r2 = 0.254m , and L = 2m . The BEM result is compared to the analytical solution in Fig. 8. For this particular configuration, the plane-wave cutoff frequency at the room temperature is 613.5Hz . It is seen that the BEM solution compares very well with the analytical solution. The minor discrepancies could be due to the collocation nature of the impedance-to-scattering matrix method and the FEM interpolation of the modes at the centroid of each element. The second test case is a rectangular bar silencer as shown in Fig. 9. This test case was first reported by Kirby et al. [27] using the pointcollocation method (hybrid analytical-FEM). With reference to Fig. 9, dx = 0.28m ，d y = 0.21m , ex = 0.06m , ey = 0.045m , and L = 0.9m . The flow resistivity for the sound absorbing material isR = 19307Rayls / m . The perforated facing sheet has porosity ϕ = 0.27, hole diameter
(15)
where
WI = 0.5S1/ ρ0 c0 n
WO = 0.5 ∑
∫ Re (ΦCi Ci+) Conj (ΦCi Ci+kCzi) dS2 /ρ0 ω
i =1
Re denotes the real part, Conj denotes the complex conjugate, n is the total number of FEM modes at the outlet, and S1 and S2 are the crosssectional areas of the inlet and outlet, respectively. 3. Test cases The first test case is a single round bar silencer housed in a rigid circular duct as shown in Fig. 7. Although this single round bar silencer does not represent a building block taken from a lattice arrangement, the analytical solution for this axisymmetric configuration can be readily obtained by using the method developed by Selamet et al. [5].
Dimensions of the Fig. 7. Configuration and dimensions of a round bar silencer.
Fig. 9. Configuration and dimensions of a rectangular bar silencer.
104
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
Fig. 11. TL prediction of the triangle module.
Fig. 10. TL comparison of the rectangular bar silencer.
dh = 0.003m , and wall thickness tw = 0.0016m . Fig. 10 compares the BEM result to the numerical solution taken from Ref. [27], along with the insertion loss (IL) measurement data from the same reference. In practice, IL is easier to measure on site, but it is not the same as the more theoretical TL. Nonetheless, it does provide a reference for comparison purposes. For large dissipative silencers, the trend of IL should be close to the trend of TL. From Fig. 10, it is seen that the proposed method compares fairly well with the point collocation method for TL up to 8000 Hz, although both numerical TL predictions are higher than the IL measurement. The third test case is a triangular module isolated from a shifted lattice arrangement of round bars shown in Fig. 4. With reference to Fig. 4, L = 2m , l1 = 0.4m , and l2 = 0.1m . Three case studies (see Table 1) with different flow resistivity and porosity combinations are used and the respective TL results are compared in Fig. 11. Increasing the porosity of the perforated facing sheet may improve the acoustic attenuation performance at higher frequencies; however, this may be at the expense of a small reduction of TL at lower frequencies. It should be noticed that at very high frequencies, the effect of flow resistivity is very little. The last test case is a comparison between the triangular module and the square module to demonstrate the advantage of the shifted lattice arrangement (Fig. 4) over the aligned lattice arrangement (Fig. 2). Fig. 12 shows the dimensions of the two different bar silencer designs. The round bar is made of polyester with flow resistivity R=16000 Rayls/m, and is covered by a 30% open perforated facing sheet. With reference to Fig. 12, l1 = 0.4m and l2 = 0.1m . The shifted lattice (triangular module) is a more compact design than the aligned lattice (square module). The TL comparison of the two designs is shown in Fig. 13, and it is seen that triangular module has better performance over 1000 Hz.
Fig. 12. Dimensions of the two different designs: (a) triangular module; (b) square module.
Fig. 13. TL comparison between the triangular and square modules.
isolated from various lattice arrangements of bar silencers are demonstrated. Although the 2D FEM is used as a post-processing tool in the modal expansion for TL computation, the proposed method is by no means a hybrid FEM-BEM or a hybrid analytical-BEM technique, as the primary analysis tool is still the same substructure BEM as in Ref. [20]. The decoupling of the main BEM solver from the post-processing TL computation is an important feature as it will allow independent modular programming and code verification. Several test cases are used to validate the proposed method. The proposed method compares very well with the analytical solution for an axisymmetric bar silencer, and with the point collocation method for a rectangular bar silencer. A triangular test case is also demonstrated and compared to the square case to show its advantage.
4. Conclusions In this paper, the 2D FEM is used to extract the cross-sectional modes numerically so that the impedance-to-scattering matrix method can be applied to any arbitrary inlet/outlet shapes. Specifically, applications of the method to the BEM modeling of a single unit Table 1 Calculation data for the triangle module. Case number
Flow resistivity (Rayls/m)
Porosity
(a) (b) (c)
16000 16000 1800
0.30 0.08 0.30
105
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
splitter silencers. J Acoust Soc Am 2005;118:2302–12. [14] Kirby R, Amott K, Williams PT, Duan W. On the acoustic performance of rectangular splitter silencers in the presence of mean flow. J Sound Vib 2014;333:6295–311. [15] Wu TW, Wan GC. Muffler performance studies using a direct mixed-body boundary element method and a three-point method for evaluating transmission loss. J Vib Acoust ASME Trans 1996;118:479–84. [16] Wu TW, Zhang P, Cheng CYR. Boundary element analysis of mufflers with an improved method for deriving the four-pole parameters. J Sound Vib 1998;217:767–79. [17] Wu TW, Cheng CYR, Zhang P. A direct mixed-body boundary element method for packed silencers. J Acoust Soc Am 2002;111:2566–72. [18] Chen JT, Chen KH. Applications of the dual integral formulation in conjunction with fast multipole method in large-scale problems for 2D exterior acoustics. Eng Anal Bound Elem 2004;28:685–709. [19] Liu Y. Fast multipole boundary element method: theory and applications in engineering. New York: Cambridge university press; 2009. [20] Lou G, Wu TW, Cheng CYR. Boundary element analysis of packed silencers with a substructuring technique. Eng Anal Bound Elem 2003;27:643–53. [21] Zhou L, Wu TW, Ruan K, Herrin DW. A reciprocal identity method for large silencer analysis. J Sound Vib 2016;364:165–76. [22] Wang P, Wu TW. Impedance-to-scattering matrix method for large silencer analysis using direct collocation. Eng Anal with Bound Elem 2016;73:191–9. [23] Åbom M. Measurement of the scattering matrix of acoustical two-ports. Mech Syst Signal Pr 1991;5:89–104. [24] Yang L, Ji ZL, Wu TW. Transmission loss prediction of silencers by using combined boundary element method and point collocation approach. Eng Anal Bound Elem 2015;61:265–73. [25] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Appl Mech Rev 1999;52:17–33. [26] Fang Z, Ji ZL. Acoustic attenuation analysis of expansion chambers with extended inlet/outlet. Noise Control Eng J 2013;61:240–9. [27] Kirby R, Williams PT, Hill J. A three dimensional investigation into the acoustic performance of dissipative splitter silencers. J Acoust Soc Am 2014;135:2727–37.
Acknowledgments The research was partially supported by the Vibro-Acoustics Consortium. L. Yang was supported by China Scholarship Council for his study at the University of Kentucky as a visiting graduate student. References [1] Munjal ML, Galaitsis AG, Vér IL. Passive silencers. In noise and vibration control engineering-principles and applications Chap. 9. In: Ver Istvan L, Beranek Leo L, editors. edited by. New Jersey: John Wiley & Sons; 2006. [2] Mechel FP. Theory of baffle-type silencers. Acustica 1990;70:93–111. [3] Mechel FP. Numerical results to the theory of baffle-type silencers. Acustica 1990;72:7–20. [4] Astley RJ, Cummings A. A finite element scheme for attenuation in ducts lined with porous material: comparison with experiment. J Sound Vib 1987;116:239–63. [5] Selamet A, Xu MB, Lee IJ. Analytical approach for sound attenuation in perforated dissipative silencers. J Acoust Soc Am 2004;115:2091–9. [6] Xu MB, Selamet A, Lee IJ, Huff NT. Sound attenuation in dissipative expansion chambers. J Sound Vib 2004;272:1125–33. [7] Venkatesham B, Munjal ML. Acoustic performance of an industrial muffler InterNoise. Portugal: Lisbon; 2010. [8] Nilsson NÅ, Söderqvist S. The bar silencer-improving attenuation by constricted two-dimensional wave propagation. Proc Inter Noise 1983;83:1–4. [9] Cummings A, Astley RJ. Finite element computation of attenuation in bar-silencers and comparison with measured data. J Sound Vib 1996;196:351–69. [10] Kirby R. A comparison between analytic and numerical methods for modeling automotive dissipative silencers with mean flow. J Sound Vib 2009;325:565–82. [11] Kirby R, Lawrie JB. A point collocation approach to modelling large dissipative silencers. J Sound Vib 2005;286:313–39. [12] Lawrie JB, Kirby R. Mode-matching without root-finding: application to a dissipative silencer. J Acoust Soc Am 2006;119:2050–61. [13] Kirby R. The influence of baffle fairings on the acoustic performance of rectangular
106
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
crossmark
Boundary element analysis of bar silencers using the scattering matrix with two-dimensional finite element modes ⁎
L. Yanga, P. Wangb, T.W. Wub, a b
School of Power and Energy Engineering, Harbin Engineering University, Harbin, Heilongjiang 150001, PR China Department of Mechanical Engineering, University of Kentucky, Lexington, KY 40506-0503, USA
A R T I C L E I N F O
A BS T RAC T
Keywords: Acoustics Mufflers Silencers Boundary element method Finite element method
Bar silencers used in industry may consist of a large array of rectangular or round bars packed in a rectangular lattice arrangement. Due to the size of the lattice, normally only a single unit that represents a building block for the lattice is isolated for analysis purposes. Even with one isolated unit, the inlet and the outlet are still quite large, and the plane-wave cutoff frequency can be very low. Therefore, higher-order modes must be considered at the inlet and outlet in order to calculate the transmission loss. This paper uses the recently developed “impedance-to-scattering matrix method” to convert the element-based impedance matrix into the mode-based scattering matrix for transmission loss calculation. Depending on the shape of the inlet and outlet, it may not always be possible to find an analytical expression of the modes needed for the modal expansion. In this paper, the two-dimensional finite element method is used to extract the eigenvalues and the eigenvectors of the inlet/ outlet cross section. The eigenvectors are then used in the modal expansion to convert the impedance matrix into the scattering matrix. Test cases include several commonly used inlet and outlet configurations, such as rectangular, circular and triangular cross sections.
1. Introduction Dissipative silencers are widely used to attenuate broadband noise generated from gas turbine engines, power plants and HVAC ducts. There are a variety of designs of dissipative silencers distributing sound absorbent in various ways. A very simple design called the lined duct design places the sound absorbing material around the circumference of the cross-sectional area. Another prevailing design is the splitter silencer or the parallel-baffle silencer, which arranges multiple sound absorbing baffles parallel to the direction of exhaust. Normally, the sound attenuation of such silencers is proportional to the perimeterarea ratio and length [1]. There have been plenty of research activities on lined-duct silencers or splitter silencers either analytically, numerically or experimentally [2–7]. In 1983, Nilsson and Söderqvist proposed the idea of bar silencers and claimed that an array of square bars made of sound absorbing materials have certain advantages over a similarly configured splitter silencer [8]. In 1996, Cummings and Astley [9] investigated the acoustical behavior of square bar silencers using the finite element method (FEM) and compared to the experimental data. To date, the FEM is still the most developed and widely used numerical technique for large dissipative silencers. Kirby and his co-
⁎
workers [10–14] used a hybrid analytical-FEM to study the acoustical performance of various large dissipative silencers. To apply the hybrid technique, the 2D FEM is first employed to extract the eigenvalues and the associated eigenvectors of an axially uniform cross section. These 2D transversal modes are then used in the modal expansion along the axial direction if the cross section remains the same. To determine the unknown amplitudes in the modal expansion, either a point collocation method or an integral-based mode matching scheme is adopted to enforce the continuity of sound pressure and particle velocity at both ends where the uniform section meets the flanges or any irregular junctions. Since the FEM is mainly used on a 2D cross section to extract the modes, the hybrid analytical-FEM is a very efficient numerical technique for silencers with a very long axially uniform section. As a viable alternative to the FEM, the boundary element method (BEM) has also been used for muffler and silencer analysis [15–17]. The direct mixed-body BEM in Refs. [15–17] can handle complex internal components, such as thin bodies, perforated tubes, and multiple bulk-reacting materials, in one single BEM domain without resorting to the multi-domain approach. It is well known that the BEM has the advantage of modeling the surface only. However, for a large silencer, even the surface mesh may still contain too many elements for
Corresponding author. E-mail address: [email protected] (T.W. Wu).
http://dx.doi.org/10.1016/j.enganabound.2016.11.001 Received 30 January 2016; Received in revised form 29 September 2016; Accepted 2 November 2016 0955-7997/ © 2016 Published by Elsevier Ltd.
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
a desktop computer to run comfortably. Normally, one way to solve large-scale problems in BEM is to use the fast multipole method (FMM) [18,19]. The FMM is especially useful for exterior radiation and scattering problems. However, the FMM has not been extended to silencer problems with complex internal components and multiple bulk-reacting materials yet. Fortunately, the interior domain of a large silencer can always be divided into several smaller substructures by using the substructure BEM [20] so that each substructure can fit within the memory limitation of a desktop computer. Each substructure produces an impedance matrix, and a so-called “impedance matrix synthesis” procedure [20] is used to combine all the substructures. Until very recently, the BEM analysis of mufflers and silencers has been limited to problems with a small inlet and a small outlet only due to its dependency on either the four-pole matrix or the anechoic termination boundary condition for transmission loss (TL) calculation [15,16]. Without a proper introduction of modal expansion into the BEM, it is difficult to model large silencers with many higher-order modes emerging at the inlet and outlet cross sections. Zhou et al. [21] recently proposed a reciprocal identity method in conjunction with the BEM impedance matrix to extract the higher-order modes at the inlet and outlet. Each reciprocal identity couples the analytical modal expansion in the inlet and outlet ducts to a BEM solution with a random boundary condition set. The BEM impedance matrix can naturally provide more than enough such solutions for the reciprocal identity coupling. The method can be regarded as an indirect postprocessing filter applied to the BEM impedance matrix to extract the higher-order modes for TL calculation. Parallel to the reciprocal identity method, a different technique called the “impedance-to-scattering matrix method” was recently developed by Wang and Wu [22] as a direct collocation approach to convert the element-based impedance matrix into the mode-based scattering matrix for TL computation. The BEM impedance matrix relates sound pressures at the inlet and outlet to the corresponding particle velocities, while the scattering matrix relates the modes at the inlet and outlet [23]. It should be noted that the term “scattering matrix” can be a little confusing and it is not the same as the T-matrix method for exterior scattering problems. The scattering matrix is a special tool designed for muffler and silencer analysis, and is equivalent to the four-pole matrix if the plane-wave assumption is still valid at the inlet and outlet. In the impedance-to-scattering matrix method [22], each sound pressure and particle velocity can be directly expanded in terms of the duct modes at the centroid of each constant boundary element. These point-wise expansions are then related by the BEM impedance matrix, and the scattering matrix can be obtained after a few matrix operations. Like the reciprocal identity method, the impedance-to-scattering matrix method can be regarded as a postprocessing filter applied to the BEM impedance matrix. The main difference between these two methods is that the reciprocal identity method is an integral-based method, while the impedance-to-scattering matrix method is a direct collocation method. Both methods rely on an analytical modal expansion at the inlet and outlet. Normally this is not an issue for a rectangular or circular inlet/outlet. However, more complicated inlet/outlet configurations will require a numerical solution of the duct modes first. In this paper, the impedance-to scattering matrix method by Wang and Wu [22] is expanded to include the use of numerical modes from the 2D FEM so that any irregular inlet/outlet configurations can be modeled. The function of the 2D FEM is to extract the eigenvalues and the associated eigenvectors (modes) of the inlet/outlet cross sections. The 2D FEM procedure is the same as the one used in the hybrid analytical-FEM [10–14] or the hybrid analytical-BEM [24], except that it is only done at the inlet/outlet cross sections as a post-processing tool to extract the duct modes. The proposed method is by no means a hybrid technique because the silencer structure itself is still modeled by the same substructure BEM as in Ref. [20] with the impedance matrix being the primary output from the BEM. Modal expansion and the
Fig. 1. Flowchart of the computation procedure.
associated 2D FEM analysis are only introduced in the end at the inlet and outlet cross sections after the resultant BEM impedance matrix is produced. Fig. 1 shows the flowchart of the proposed computation procedure. Figs. 2–4 show three typical lattice arrangements of bar silencers. Usually a small building block is isolated from the lattice for analysis purposes. The building block in Figs. 2 and 3 is simply a round or rectangular bar housed in a rectangular duct. The duct walls are assumed rigid due to symmetry. The building block in Fig. 4 actually begins with a hexagon duct. Due to rotational symmetry, the hexagon eventually reduces to just a triangle. Although rectangular and circular
Fig. 2. A rectangular module isolated from an aligned lattice arrangement of round bars.
101
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
subdomain integral equations are summed together, while Eq. (2) is the normal-derivative integral equation derived from Eq. (1).
⎞ ∂G 0 (rp, rQ ) + + jk 0 z 0 vn (rQ ) G0 (rp,rQ ) ⎟ dS + ∫ (p − p− ) dS ∂nQ T ⎠ ∼ ∼ ⎞ ⎛ ∂G (rp, rQ ) ∂G (rp, rQ ) + ∼ ∼ + ∫ ⎜p (rQ ) ∂n + jk z∼vn (rQ ) G (rp,rQ ) ⎟ dS + ∫ (p − p− ) dS ∂nQ B⎝ BTB Q ⎠ ∼ ∼ ⎡ ⎛ ∂G 0 (rp, rQ) ∂G (rp, rQ ) ⎞ ∂G (rp, rQ ) ∼ + ∫ ⎢p0 ⎜ ∂n − ∂n + jk 0 z 0 vn (rQ ) G0 (rp, rQ ) ⎟ − z 0 ζ vn (rQ ) ∂n IP ⎣ ⎝ Q Q Q ⎠ ⎛
∫R ⎜p (rQ ) ⎝
∂G 0 (rp, rQ ) ∂nQ
⎛ ∂G 0 (rp, rQ) ∼ ∂G∼(rp, rQ) ⎞ ∼ ∼ − jk z∼vn (rQ ) G (rp, rQ )] dS + ∫ ⎜p0 ∂n − p ∂n ⎟ dS ATB ⎝ Q Q ⎠ ⎧ p (rp ), P ∈ Ω. ⎪ ⎪ 0.5p (rp ), P ∈ R + B. ⎪ + − = ⎨ 0.5[p (rp ) + p (rp )], P ∈ T + BTB. ⎪ ∼ ⎪ p0 (rp ) + 0.5z 0 k 0 ζ vn (rp ), P ∈ IP. ⎪ 0.5[p (rp ) + ∼ p (rp )], P ∈ ATB. ⎩ 0
Fig. 3. A rectangular module isolated from an aligned lattice arrangement of square bars.
⎛
∂ 2G 0 (rp, rQ )
⎝
∂nQ ∂np
∫R ⎜p (rQ )
∂G 0 (rp,rQ ) ⎞
⎟ dS + ∫ T ⎠ ∼ ⎞ ∂G (rp,rQ ) ∼ + jk z∼vn (rQ ) ∂n ⎟ dS p ⎠
+ jk 0 z 0 vn (rQ )
∼ ⎛ ∂ 2G (rp, rQ ) + ∫ ⎜p (rQ ) ∂n ∂n Q p B⎝ ⎡ ⎛ ∂ 2G 0 (rp, rQ) + ∫ ⎢p0 ⎜ ∂n ∂n − Q p IP ⎣ ⎝
∂ 2G 0 (rp, rQ )
∂np
∂nQ ∂np
(1)
(p+ − p− ) dS
∼ ∂ 2G (rp, rQ ) ⎞ ∂nQ ∂np
⎟ ⎠
∼ ∼ ∂ 2G (rp, rQ ) ∂G 0 (rp, rQ ) ∂G (rp, rQ ) ⎤ ∼ ∼ − z 0 ζ vn (rQ ) ∂n ∂n + jk 0 z 0 vn (rQ ) ∂n −jk z∼vn (rQ ) ∂n ⎥ dS Q p p p ⎦
⎛ ∂ 2G 0 (rp, rQ) ∼ ∂ 2G∼(rp, rQ) ⎞ + ∫ ⎜p0 ∂n ∂n − p ∂n ∂n ⎟ dS + ∫ BTB ATB ⎝ Q p Q p ⎠ ⎧ 0, P ∈ T + ATB + BTB. =⎨ ∼∼ ⎩−0.5(jk 0 z 0 + jk z ) vn (rp ), P ∈ IP.
∼ ∂ 2G (rp, rQ ) ∂nQ ∂np
(p+ − p− ) dS
⎪ ⎪
(2) ∼ ∼ − jkr PQ where G0 = and G = e /4πrPQ are the free-space Green's functions in the air and the sound-absorbing material, respectively, P and Q are the collocation point and the integration point, respectively, nP and nQ are the corresponding normal vectors, rPQ = rP − rQ , z 0 and z∼ represent the characteristic impedances of air ∼ and the sound-absorbing material, k 0 and k the wavenumbers of ∼ corresponding media, respectively, ζ the transfer impedance of any perforated interfaces. The symbols R, T, B, IP, ATB and BTB represent “regular”, “thin”, “bulk-reacting”, “interface with a perforated facing sheet”, “air-thin-bulk reacting” and “bulk reacting-thin-bulk reacting” surfaces as defined in Fig. 5. p denotes sound pressure (on IP or ATB, ∼ are the sound pressure variables in the air and the soundp0 and p absorbing material, respectively) and vn is the normal particle velocity. The unit normal vector on either R or B is pointing into the interior acoustic domain. On T and BTB, the unit normal vector can be pointing into either side of the thin surface. The side into which the normal is pointing is called the positive side while the opposite side is called the negative side. The corresponding sound pressure values are denoted by p+ and p−, respectively. The unit normal vector on IP or ATB is pointing
e−jk 0 rPQ /4πrPQ
Fig. 4. A triangular module isolated from a shifted lattice arrangement of round bars.
duct modes are well known analytically, triangular or hexagon duct modes may not be easy to find. In some single-silencer applications, the inlet and/or outlet may also be an ellipse, trapezoid, or other irregular shapes. In this paper, the 2D FEM is used to extract the duct modes at the inlet and outlet, regardless of the shape. The method is first verified by comparing to an analytical solution for a single round bar silencer housed in a circular duct. The result of a rectangular bar-silencer test case is validated by comparing to an existing hybrid analytical-FEM solution available in literature. Then a triangular test case follows and is compared to the rectangular design.
2. Theory 2.1. Direct mixed-body BEM The concept of the direct mixed-body BEM [15–17,20] actually begins with the conventional multi-domain BEM by subdividing the acoustic domain into several well-defined homogeneous subdomains. The Helmholtz integral equation may be written for each individual subdomain. All the integral equations are then summed to create a single integral equation. The normal-derivative hypersingular integral equation is used to provide an additional equation at any interfaces that have two unknown variables. A review of hypersingular integrals can be found in Chen and Hong [25]. The direct mixed-body BEM is ideally suited to mufflers and silencers with complex internal components. The integral equations needed in this paper are listed in Eq. (1) and Eq. (2) and the complete edition can be found in Ref. [20]. Eq. (1) is the conventional Helmholtz integral equation except that all the
Fig. 5. Types of surfaces defined in the direct mixed-body BEM.
102
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
circular or rectangular inlet/outlet, the eigenfunctions can be obtained analytically. For irregular inlet/outlet shapes, the eigenvalues and the corresponding eigenfunctions must be obtained numerically. In the paper, the 2D FEM is used to calculate the eigenvalues and the associated eigenfunctions of the inlet/outlet cross sections. The governing differential equation for the transversal sound pressure in a 2D cross section is [26] 2 ∇2xy pxy + kxy pxy = 0
(6)
where pxy is the transversal sound pressure, and kxy is the transversal wavenumber that must satisfy the following equation: 2 kxy + kz2 = k 02
In FEM, the transversal sound pressure at any point is expressed in terms of the shape functions
Fig. 6. Substructure technique.
pxy = N Tp
into the air side.
(8)
where the matrix N consists of the column vectors of the global shape functions, and p is the column vector of the nodal values of transversal sound pressure. Applying the rigid-wall boundary condition and the Galerkin procedure [26] to Eq. (6) yields the following eigen-equation:
2.2. Impedance matrix and substructure BEM For large silencers, the substructure BEM [20] is a very useful technique as it can be difficult to analyze a large silencer in one single BEM model on a desktop computer. A demonstrated in Fig. 6, a large silencer can be divided into several smaller substructures so that each substructure can fit within the memory limitation of a desktop computer. Specifically, if a silencer has an axially uniform section, the impedance matrix of a small template taken from the uniform section can be repeatedly used without having to model the long uniform section. The impedance matrix of each substructure relates sound pressures at the inlet and outlet of the substructure to the corresponding particle velocities. Multiple computers may be used simultaneously for different substructures as each substructure analysis is independent of the others. All the substructure impedance matrices are then combined to form the resultant impedance matrix for the whole silencer. The detailed procedure is described in Ref. [20]. The end result is the resultant impedance matrix for the large silencer:
⎡ p I ⎤ ⎡ Z11 Z12 ⎤ ⎡ vI ⎤ ⎥⎢ ⎥ ⎢p ⎥ = ⎢ ⎣ O ⎦ ⎣ Z21 Z22 ⎦ ⎣ vO ⎦
(7)
2 (K − kxy M) p = 0,
(9)
(∇N)e (∇N)Te dSe Se
where K = ∑ ∫
(N)e (N)Te dSe Se
and M = ∑ ∫
are the
stiffness matrix and the mass matrix, respectively, subscript "e" denotes the element. If the number of nodes is n, n eigenvalues (kxy )i (1 ≤ i ≤ n ) and the associated eigenvectors (Φ)i with length of n may be obtained by solving Eq. (9). Assemble all the eigenvalues and eigenvectors in the vector k and the matrix Φ , respectively. Let A denote the wave amplitudes at the inlet and C the wave amplitudes at the outlet. Apply Eqs. (4) and (5) to the inlet and outlet individually, and sound pressures and particle velocities at the inlet and the outlet may be written in terms of their corresponding modes as
(3)
where the vectors p and v represent the sound pressures and particle velocities, and the subscripts I and O denote the inlet and the outlet, respectively.
ΦAn ⎤ ⎡ A+ ⎤ ⎡ p I ⎤ ⎡ ΦAn n ⎢ v ⎥ = ⎢ kAzn Φ − kAzn Φ ⎥ ⎢ − ⎥ ⎣ I ⎦ ⎢⎣ ρ ω An ⎥ ρ0 ω An ⎦ ⎣ An ⎦ 0
(10)
ΦCn ⎤ ⎡ C+ ⎤ ⎡ pO ⎤ ⎡ ΦCn ⎥⎢ n ⎥ k ⎢ v ⎥ = ⎢ kCzn Φ ΦCn ⎥ ⎣ C−n ⎦ ⎣ O ⎦ ⎢⎣ ρ ω Cn − ρCzn ⎦ 0 0ω
(11)
where kAzn and k Czn contain all the eigenvalues, and ΦAn and ΦCn contain the corresponding eigenvectors. From Eqs. (3), (10), and (11), the following relationship is obtained
2.3. Impedance-to-scattering matrix method with FEM modes
⎤ ⎧ C+n ⎫ ⎡ MPI−Z11MVI ⎤ ⎧ A+n ⎫ ⎡ Z12 MVO ⎥⎨ ⎬ ⎢ ⎥ ⎨ −⎬ = ⎢ ⎣ Z21MVI ⎦ ⎩ An ⎭ ⎣ MPO−Z22 MVO ⎦ ⎩ C−n ⎭
The impedance matrix in Eq. (3) can be converted into the scattering matrix by expressing each sound pressure and particle velocity at the inlet and outlet in terms of the duct modes [22]. In this paper, the duct modes are obtained by the 2D FEM so that the impedance-to-scattering matrix method is not limited to certain regular shapes only. At the inlet and outlet, sound pressure and particle velocity at any point nC can be expressed in a modal expansion:
(12)
where
MPI = [ ΦAn ΦAn ] MPO = [ ΦCn ΦCn ]
⎡k ΦAn − MVI = ⎢ ρAzn ⎣ 0ω
⎤ kAzn Φ ρ0 ω An ⎥ ⎦
M
pnC =
∑ ΦnC m (Fm+e−jkzm z + Fm−e jkzm z ) m =1
⎡k ΦCn − MVO = ⎢ ρCzn ⎣ 0ω
(4)
and
vnC =
1 ρ0 ω Fm+
Although Eq. (12) may look like a “transfer scattering matrix” that relates the wave amplitudes at the inlet to the wave amplitudes at the outlet, matrix inverse should not be performed at this stage because Eq. (12) does not represent a well-posed boundary value problem. A wellposed boundary value problem should have a known condition at the inlet and another known condition at the outlet. For example, to find the TL, a certain incident wave should be given at the inlet and the outlet is assumed anechoic. Rearrange Eq. (12) to get
M
∑ ΦnC m kzm (Fm+e−jkzm z − Fm−e jkzm z ) m =1
⎤ k Czn ΦCn ⎥ ρ0 ω ⎦
(5)
Fm−
where and are the modal amplitudes corresponding to the acoustic waves of order m travelling in the positive and negative z directions, respectively, kz represents the wavenumber in the axial direction, ΦnC m denotes the eigenfunction value at the point nC . For a 103
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
⎧ C+ ⎫ ⎧ A+ ⎫ [SO] ⎨ −n ⎬ = [SI] ⎨ −n ⎬ ⎩ An ⎭ ⎩ Cn ⎭
(13)
Since there are always more boundary elements at the inlet and outlet than the number of modes, a least-square matrix inverse may be performed on Eq. (13) to obtain the scattering matrix S :
⎧ C+n ⎫ ⎧ A+ ⎫ ⎨ − ⎬ = [S] ⎨ −n ⎬ ⎩ An ⎭ ⎩ Cn ⎭
(14)
2.4. Transmission loss To calculate the transmission loss, a single incident plane wave is prescribed at the inlet, and the outlet is assumed anechoic. In other words, let A+n = [1, 0, 0, ... ,0]T and C−n = [0, 0, 0, ... ,0]T . The transmitted wave amplitudes C+n and the reflected wave amplitudes A−n can be obtained from the scattering matrix relationship, Eq. (14). The transmission loss is defined as the difference of sound power in dB between the single incident plane wave and the transmitted waves. That is,
TL = 10 log10 (WI / WO )
Fig. 8. TL comparison of the round bar silencer.
It is noted that a similar test case was also used in Refs. [21,22] to validate the reciprocal identity method and the impedance-to-scattering matrix method, respectively. The purpose of repeating this particular test case here is to validate the proposed scattering matrix conversion method using the 2D FEM modes. The round bar is made of polyester with flow resistivity R=16000 Rayls/m, and is covered by a 30% open perforated facing sheet to protect the material from being blown away by the exhaust gas. Both ends of the bar are covered by a steel plate. With reference to Fig. 6, the dimensions used in this test case are r1 = 0.342m , r2 = 0.254m , and L = 2m . The BEM result is compared to the analytical solution in Fig. 8. For this particular configuration, the plane-wave cutoff frequency at the room temperature is 613.5Hz . It is seen that the BEM solution compares very well with the analytical solution. The minor discrepancies could be due to the collocation nature of the impedance-to-scattering matrix method and the FEM interpolation of the modes at the centroid of each element. The second test case is a rectangular bar silencer as shown in Fig. 9. This test case was first reported by Kirby et al. [27] using the pointcollocation method (hybrid analytical-FEM). With reference to Fig. 9, dx = 0.28m ，d y = 0.21m , ex = 0.06m , ey = 0.045m , and L = 0.9m . The flow resistivity for the sound absorbing material isR = 19307Rayls / m . The perforated facing sheet has porosity ϕ = 0.27, hole diameter
(15)
where
WI = 0.5S1/ ρ0 c0 n
WO = 0.5 ∑
∫ Re (ΦCi Ci+) Conj (ΦCi Ci+kCzi) dS2 /ρ0 ω
i =1
Re denotes the real part, Conj denotes the complex conjugate, n is the total number of FEM modes at the outlet, and S1 and S2 are the crosssectional areas of the inlet and outlet, respectively. 3. Test cases The first test case is a single round bar silencer housed in a rigid circular duct as shown in Fig. 7. Although this single round bar silencer does not represent a building block taken from a lattice arrangement, the analytical solution for this axisymmetric configuration can be readily obtained by using the method developed by Selamet et al. [5].
Dimensions of the Fig. 7. Configuration and dimensions of a round bar silencer.
Fig. 9. Configuration and dimensions of a rectangular bar silencer.
104
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
Fig. 11. TL prediction of the triangle module.
Fig. 10. TL comparison of the rectangular bar silencer.
dh = 0.003m , and wall thickness tw = 0.0016m . Fig. 10 compares the BEM result to the numerical solution taken from Ref. [27], along with the insertion loss (IL) measurement data from the same reference. In practice, IL is easier to measure on site, but it is not the same as the more theoretical TL. Nonetheless, it does provide a reference for comparison purposes. For large dissipative silencers, the trend of IL should be close to the trend of TL. From Fig. 10, it is seen that the proposed method compares fairly well with the point collocation method for TL up to 8000 Hz, although both numerical TL predictions are higher than the IL measurement. The third test case is a triangular module isolated from a shifted lattice arrangement of round bars shown in Fig. 4. With reference to Fig. 4, L = 2m , l1 = 0.4m , and l2 = 0.1m . Three case studies (see Table 1) with different flow resistivity and porosity combinations are used and the respective TL results are compared in Fig. 11. Increasing the porosity of the perforated facing sheet may improve the acoustic attenuation performance at higher frequencies; however, this may be at the expense of a small reduction of TL at lower frequencies. It should be noticed that at very high frequencies, the effect of flow resistivity is very little. The last test case is a comparison between the triangular module and the square module to demonstrate the advantage of the shifted lattice arrangement (Fig. 4) over the aligned lattice arrangement (Fig. 2). Fig. 12 shows the dimensions of the two different bar silencer designs. The round bar is made of polyester with flow resistivity R=16000 Rayls/m, and is covered by a 30% open perforated facing sheet. With reference to Fig. 12, l1 = 0.4m and l2 = 0.1m . The shifted lattice (triangular module) is a more compact design than the aligned lattice (square module). The TL comparison of the two designs is shown in Fig. 13, and it is seen that triangular module has better performance over 1000 Hz.
Fig. 12. Dimensions of the two different designs: (a) triangular module; (b) square module.
Fig. 13. TL comparison between the triangular and square modules.
isolated from various lattice arrangements of bar silencers are demonstrated. Although the 2D FEM is used as a post-processing tool in the modal expansion for TL computation, the proposed method is by no means a hybrid FEM-BEM or a hybrid analytical-BEM technique, as the primary analysis tool is still the same substructure BEM as in Ref. [20]. The decoupling of the main BEM solver from the post-processing TL computation is an important feature as it will allow independent modular programming and code verification. Several test cases are used to validate the proposed method. The proposed method compares very well with the analytical solution for an axisymmetric bar silencer, and with the point collocation method for a rectangular bar silencer. A triangular test case is also demonstrated and compared to the square case to show its advantage.
4. Conclusions In this paper, the 2D FEM is used to extract the cross-sectional modes numerically so that the impedance-to-scattering matrix method can be applied to any arbitrary inlet/outlet shapes. Specifically, applications of the method to the BEM modeling of a single unit Table 1 Calculation data for the triangle module. Case number
Flow resistivity (Rayls/m)
Porosity
(a) (b) (c)
16000 16000 1800
0.30 0.08 0.30
105
Engineering Analysis with Boundary Elements 74 (2017) 100–106
L. Yang et al.
splitter silencers. J Acoust Soc Am 2005;118:2302–12. [14] Kirby R, Amott K, Williams PT, Duan W. On the acoustic performance of rectangular splitter silencers in the presence of mean flow. J Sound Vib 2014;333:6295–311. [15] Wu TW, Wan GC. Muffler performance studies using a direct mixed-body boundary element method and a three-point method for evaluating transmission loss. J Vib Acoust ASME Trans 1996;118:479–84. [16] Wu TW, Zhang P, Cheng CYR. Boundary element analysis of mufflers with an improved method for deriving the four-pole parameters. J Sound Vib 1998;217:767–79. [17] Wu TW, Cheng CYR, Zhang P. A direct mixed-body boundary element method for packed silencers. J Acoust Soc Am 2002;111:2566–72. [18] Chen JT, Chen KH. Applications of the dual integral formulation in conjunction with fast multipole method in large-scale problems for 2D exterior acoustics. Eng Anal Bound Elem 2004;28:685–709. [19] Liu Y. Fast multipole boundary element method: theory and applications in engineering. New York: Cambridge university press; 2009. [20] Lou G, Wu TW, Cheng CYR. Boundary element analysis of packed silencers with a substructuring technique. Eng Anal Bound Elem 2003;27:643–53. [21] Zhou L, Wu TW, Ruan K, Herrin DW. A reciprocal identity method for large silencer analysis. J Sound Vib 2016;364:165–76. [22] Wang P, Wu TW. Impedance-to-scattering matrix method for large silencer analysis using direct collocation. Eng Anal with Bound Elem 2016;73:191–9. [23] Åbom M. Measurement of the scattering matrix of acoustical two-ports. Mech Syst Signal Pr 1991;5:89–104. [24] Yang L, Ji ZL, Wu TW. Transmission loss prediction of silencers by using combined boundary element method and point collocation approach. Eng Anal Bound Elem 2015;61:265–73. [25] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Appl Mech Rev 1999;52:17–33. [26] Fang Z, Ji ZL. Acoustic attenuation analysis of expansion chambers with extended inlet/outlet. Noise Control Eng J 2013;61:240–9. [27] Kirby R, Williams PT, Hill J. A three dimensional investigation into the acoustic performance of dissipative splitter silencers. J Acoust Soc Am 2014;135:2727–37.
Acknowledgments The research was partially supported by the Vibro-Acoustics Consortium. L. Yang was supported by China Scholarship Council for his study at the University of Kentucky as a visiting graduate student. References [1] Munjal ML, Galaitsis AG, Vér IL. Passive silencers. In noise and vibration control engineering-principles and applications Chap. 9. In: Ver Istvan L, Beranek Leo L, editors. edited by. New Jersey: John Wiley & Sons; 2006. [2] Mechel FP. Theory of baffle-type silencers. Acustica 1990;70:93–111. [3] Mechel FP. Numerical results to the theory of baffle-type silencers. Acustica 1990;72:7–20. [4] Astley RJ, Cummings A. A finite element scheme for attenuation in ducts lined with porous material: comparison with experiment. J Sound Vib 1987;116:239–63. [5] Selamet A, Xu MB, Lee IJ. Analytical approach for sound attenuation in perforated dissipative silencers. J Acoust Soc Am 2004;115:2091–9. [6] Xu MB, Selamet A, Lee IJ, Huff NT. Sound attenuation in dissipative expansion chambers. J Sound Vib 2004;272:1125–33. [7] Venkatesham B, Munjal ML. Acoustic performance of an industrial muffler InterNoise. Portugal: Lisbon; 2010. [8] Nilsson NÅ, Söderqvist S. The bar silencer-improving attenuation by constricted two-dimensional wave propagation. Proc Inter Noise 1983;83:1–4. [9] Cummings A, Astley RJ. Finite element computation of attenuation in bar-silencers and comparison with measured data. J Sound Vib 1996;196:351–69. [10] Kirby R. A comparison between analytic and numerical methods for modeling automotive dissipative silencers with mean flow. J Sound Vib 2009;325:565–82. [11] Kirby R, Lawrie JB. A point collocation approach to modelling large dissipative silencers. J Sound Vib 2005;286:313–39. [12] Lawrie JB, Kirby R. Mode-matching without root-finding: application to a dissipative silencer. J Acoust Soc Am 2006;119:2050–61. [13] Kirby R. The influence of baffle fairings on the acoustic performance of rectangular
106