- Email: [email protected]
A study of sound transmission in curved duct bends by the Galerkin method
Journal of Sound and Vibration (1976) 45(l),
91-104
A STUDY OF SOUND TRANSMISSION DUCT BENDS BY THE GALERKIN
IN CURVED METHOD
C. K. W. TAM Department of Mathematics, Florida State University, Tallahassee, Florida 32306, U.S.A.
(Received 19 June 1975, and in revisedform
8 August 1975)
A computational procedure based on the Galerkin method is developed to study the sound transmission and reflection characteristics of circularly curved bends of rectangular ducts. This procedure is adaptable to routine computer calculation. An important component of the computer program consists of a QR iterative algorithm which solves the matrix eigenvalue problem generated by the Galerkin method. The present procedure produced very satisfactory convergent results when applied even to cases where the incident wave has high frequency and high wave mode numbers. Some properties of the Galerkin solution are investigated. Its relationship to the classical solution by the method of separation of variables is discussed. 1. INTRODUCTION
This paper is concerned with the computational aspect of the transmission of sound through circularly curved bends of a duct system with rectangular cross-section. Curved bends are unavoidable features of most practical duct systems. An efficient procedure to compute their acoustical characteristics should, therefore, be valuable to engineers and architects. The propagation of sound waves in curved ducts is characterized by wave patterns which differ in many important respects from those in straight ducts. In a duct system consisting of straight ducts linked together by curved bends the continuous reflection of sound waves from the curved boundary of a bend causes a partial reflection of any incident wave from the adjoining straight ducts. This reduces the sound transmission capability of the duct system as compared with that of a single straight duct. The physics of the problem is rather simple and can easily be understood. The main problem from a practical standpoint is computational difficulty. For a circular bend of a rectangular duct it was recognized by Krasnushkin [l] and others, some years ago, that the equations of motion were separable with respect to a suitably chosen cylindrical co-ordinate system. However, for a bend of finite length this solution must be joined to the incident, transmitted and reflected wave solutions of the adjoining straight ducts. To avoid these problems Krasnushkin and later Seveshnikov [2], Grigor’yan [3] and others restricted themselves to the consideration of infinitely long bends. This is, of course, an idealization and it corresponds approximately to that of a helical wave guide. Physically, their findings could therefore provide only qualitative information for application to practical duct systems. Seveshnikov [2] generalized the problem considered by Krasnushkin [I]. He investigated theoretically acoustic wave propagation in arbitrarily deformed wave guides. On starting from the acoustic wave equation written for a general polar co-ordinate system he studied solutions which would correspond to a helical elbow bend of a circular section pipe. Unfortunately, the validity of his results is restricted to ducts with bends of very large radius of 91
92
C. K. W. TAM
curvature to diameter ratios. Osborne [4] investigated theoretically the separable solution of sound wave propagation in rectangular section curved ducts. The radial equation leads to Bessel functions of non-integer order which are exceedingly difficult to compute. He limited his studies to cases of very large inner and outer bend radii and very small inner radius. Recently works on the subject of sound transmission in curved rectangular ducts were published by Rostafinski [5, 61, Cummings [7] and Myers and Mungur [8]. Their approach follows very much the method of separation of variables used by previous investigators. However, Rostafinski and Cummings did make an effort to join the solution in the curved part of the duct to the propagating wave solutions in the straight portion of the duct system. Again the major problem of having to deal with Bessel functions of non-integer order forced Rostafinski to consider only very low frequency sound and Cummings to study only the propagation of the lowest (plane) wave mode. In his more recent publication Rostafinski [6] tried to avoid the problem of computing Bessel functions of non-integer order by treating the problem in an inverse manner. In this way he was able to produce some numerical results including velocity distributions across the duct. Unfortunately, this method does not guarantee that the full set of eigenfunctions are being used. Thus, it is possible that the computed results are incomplete. Moreover, this method, even if it is free of the above difficulty, is very cumbersome to use in a practical situation because it does not deal with the given physical problem directly. The primary objective of this investigation is to develop a direct, reliable computational procedure which can be used to determine the sound transmission and reflection characteristics of circular bends of rectangular ducts. The procedure is to be adaptable to routine computer calculations so that it can provide useful information for design or be used as a guide for other engineering purposes where sound transmission is a factor to be taken into consideration. For this purpose, the method of separation of variables which has been used by all previous investigators is found to be inadequate. In this paper the method of Galerkin modified to suit the present problem is adopted instead. The Galerkin method has been used by many authors to provide solutions to various engineering problems. A description of this method can be found in books (e.g., that of Finlayson [9]) and will not be repeated here. In duct acoustics, this method has recently been employed successfully by Eversman, Cook and Beckemeyer [lo] to study the propagation of sound in ducts with variable cross-section. For the problem under consideration, the Galerkin method, as will be shown in section 2, reduces the wave equation to a large system of linear ordinary differential equations with constant coefficients. It is well known that the solution of these equations can be constructed in terms of the eigenvectors and eigenvalues of the coefficient matrix. Thus, the success of the Galerkin method hinges on how efficiently the matrix eigenvalues problem can be solved numerically. Fortunately, there is available a stable iterative QR algorithm (see the book by Wilkinson [I 11) which is particularly effective in handling this aspect of the computation. The combined Galerkin method-QR algorithm has now been programmed into a special subroutine at the Florida State University Computer Center. Through direct numerical experiments using the Center’s CDC 6500 Computer, it is found that the procedure is well suited for the goals stated. The cases tested included high frequency, high order wave modes incident on a sharp circular bend. These tests are believed to cover not only most of the practical cases but also a wide range of parameters for which the method of separation of variables used previously failed to produce fruitful numerical results. Some properties of the solution obtained by the Galerkin method will be presented in section 3. These properties are found to be very helpful in interpreting the solution physically and mathematically. Further, they provide certain information which has proved to be valuable in actual numerical computations. Some experience gained in applying the present computation procedure will be discussed at the end of this paper.
SOUND TRANSMISSION IN DUCT BENDS
93
2. FORMULATION
Consider the transmission of sound through a circularly curved section of a duct system as shown in Figure 1. Here it is assumed that the cross-section of the duct is the same on both
I
Figure 1. Curved duct geometry.
sides of the bend. Without loss of generality the problem is conceived to consist of an incident wave with mode numbers (j,k) propagating in the positive z-direction along region (1) of the duct system. On reaching the curved portion of the duct, namely region (2), part of the wave would be reflected. The reflected wave propagates in the negative z-direction away from the x-y plane. The other part of the incident wave is transmitted into region (3) and radiates away from the c-u plane. Let p, U, u and w denote the pressure and the velocity components in the x, y and z, or 5, y and q (see Figure 2), directions, respectively. In terms of
7
Figure 2. Co-ordinate systems used.
dimensionless variables with length scale b (width of duct in the x-direction), velocity scale c (the speed of sound) time scale b/c and pressure scale pc2 (p is the gas density), the incident, reflected and transmitted waves in regions (1) and (3) can be written as follows : the incident wave (j, k)t mode : P (i)=p~)cos[lc(j2.P = i{Tc(jvti) = i{$(k
I)x]cos[nj?(k-
l)/~}p:~)sin[7c(j- l)/o}pj”
w(*) = (~~/o)p~~)cos[dj-
cos[+-
l)y]e’(KJL-W’), 1) x] cos[$(k l)x]sin[$(k
- l)y] e’(KJ’-“f), - l)y]e’(KJz-Or,r),
1)x] cos[nfi(k - I)y]ei(“Jz-‘Ot),
(1)
t Throughout this paper the wave modes are numbered in such a way that the plane wave mode is represented byj=l,k=l.
94
C. K. W. TAM
where ~~ = [co2 - n'(j- 1)”- rc2(k- l)*/?*]i/*, p,w is the amplitude of the incident wave and B = b/h is the aspect ratio of the duct; the reflected waves: P (r) =f&p;’
cos[7r(m -
1) X]cos[rcfi(k - l)y] e-i(Kmr+~t),
u(‘) =mzl i{rc(m - l)/o}pg) sin[rr(m - 1) X] cos[r$(k - l)y] e-i(Kmz+wt), n(r) = $i i{@(k - 1)/o] P,(‘) cos[lr(m
-
1) x] sin[ nj3(k - 1) y] e--i(Kmz+of),
WC’=) $,-(Km/cO)P(m’)cos[rc(rr?- 1) x] cos[rcj?(k - l)y] e-i(K~z+of),
(2)
where if w2 > n’(m - l)* + x2 /I*& - l)z,
K, = [W’- 7?()72- 1)2 - 7r2jY2(k- 1)2]1/2; ( i[n*(m - l)* + 7r2j12(k - 1)2 - ~0~1~‘~;
if o2 < 7r*(m- 1)2 + 7c2j12(k - 1)2;
the transmitted waves: P (f) =~lpg)
cos[7c (m - 1) [] cos[@(k - l)y] e’(xmn--of),
u(‘) =zii(a(m
- 1)/w} py’ sin[7r(m - 1) c] cos[n/?(k - 1) y] ei(emq-a)r),
- 1) 51sin[ a/?+ - 1) y] ei(Km*-aof), u(‘) =mzl i{rr(k - 1) /3/o}ps’ cos[ ~~(112
w(‘)=~l(K”/4P:)
cos[ 7t(m - 1) [] cos[ aj?(k - 1) y] ei(~mU-or).
(3)
Equations (I), (2) and (3) are, of course, solutions of the equations of motion,
ap _67+ r*v=o,
$&VP.
(4)
They satisfy the rigid wall boundary conditions at the surface of the duct. In addition, the radiation condition has been invoked in selecting the form of the reflected and transmitted waves. The amplitudes (complex in general) p$) and pg’ are unknown at this stage. The determination of their relationship to pj*) is the central problem of this work. In the curved part of the duct system i.e., region (2) of Figure 1 and Figure 2, the equations of motion in cylindrical co-ordinates are
ap au i aw u -+-+--+-+-~~, at ar r ae r ;+g=o,
au ay
au ap -+--0, at ay
fY+!&o.
(5)
A solution of equations (5) will now be constructed by the Gale&in method. According to this method the dependent variables are to be represented by a complete set of basis functions with arbitrary amplitudes. Here, the following basis functions will be used : p = J,
p,(O) cos[w(m - 1) (r - R)] cos[nfi (k - 1) y] e-‘“‘,
u = mz, i{lr(m - l)/o}p,(8)
sin[rc(m - l)(r - R)] cos[r$(k - 1) y] e-i”,
SOUND TRANSMISSION
v
=mzl i(7c(k-
1) /?/o}pJe)
95
IN DUCT BENDS
cos[rr(m - l)(r - R)] sin[7$(k - l)y] e-‘“‘,
w = $, (i/w) w,(8) cos[rr(m - 1) (r - R)] cos[r$(k - 1) y] e-rat.
(6)
Here R Q r < R + I, 0 < y < l/p, R being the inner radius of curvature of the curved duct. In principle, there is very little restriction as to what type of function could be used as basis functions so long as the functions form a complete set preferably satisfying certain boundary conditions. However, it is well known that for a given problem some basis functions work better than others. For this reason, it is always advantageous to use basis functions that are physically or mathematically natural to the problem. The basis functions of equations (6) are chosen for two reasons. First, they are simple and are closely related to the duct modes of a straight duct. Secondly, at f3= 0 and 0 = eE these functions match automatically with the incident, reflected and transmitted wave modes. This latter property will facilitate the process of joining these solutions together at these two planes. The two sets of functions Pm(B)and w,(e) are as yet arbitrary. They are to be chosen such that the continuity and &momentum equations of equations (5) are satisfied approximately. This is done by first substituting equations (6) into the first and fourth of equations (5) to obtain the residuals. Then the residuals are forced to become zero by requiring them to be orthogonal to a complete set of functions. For this purpose the following complete set of functions will be used : (7) r cos[n(n - 1) (r - R)], n=l,2,3,... The residuals are multiplied by each of the functions (7) and set equal to zero after integration over r from R to R + 1. In this way two infinite sets of equations for pm(e) and w,,,(e)are found. In practice these systems of equations and also the summations in equations (2), (3) and (6) are all truncated at a large integer value : say, M = n = N. The resulting finite system of equations and unknowns are hence rendered solvable. On carrying out the above steps it is easy to find that&@ and w,(e) are given by the solution of the following systems of ordinary differential equations with constant coefficients. dw,/de = 5 A,,P,,
dp,W’ = - 5 B,,, w,,
(n
=
1,2,3,...N)
where ;
&m =
[(-l).+,-
l]
[(m - 1)2 + (n - 1)2] [(m - I)2 - (n - l)‘]’ ’
n#m
_
I
n=rn
;+R, [o2
A,,,, =
(8)
m=l
m=l
_
x2
B2(k
_
1)21
B”,
_4(m- II2(n- II2rel)m+n- 11) [(m - 1)2 - (n - 1)2]2
[02-~2j?2(k-1)2-~2(m-l)2](1/2+R),
n#m n=m.
To facilitate the solution of equations (8) the following matrix notation will be introduced A = matrix [A,,], B = matrix [B,,,], P = column vector with elements p#), w = column vector with elements w,(e), p(‘) = column vector with elements p’,), p@)= column vector with elements pz), e, = column vector whosejth element is unity and all the other elements are zero.
96
c. K. w. TAM
It is important to note that the matrices A and B are real and symmetric. In matrix notation equations (8) can be rewritten as dP/dQ = -Bw.
dw/d0 = AP,
(9)
On eliminating w from equations (9), a single equation for P is found: d2 P/de2 = -BAP.
(10)
A general solution of equation (10) can easily be constructed in terms of the eigenvalues and eigenvectors of the matrix -BA. Suppose A:, A:. . . A: are the eigenvalues and CQ,a,, . . . aN are the corresponding eigenvectors of -BA: i.e., -BAa, = 1: a,.
(11)
The following matrices then can be formed : X = [al, a2,. . . aN], a matrix whose ith column is ai, l/A1
0
0
0
l/l2
0
0
l/A,
A=0
. . 0
r
eAle 0 0. . . . 0
-
0
:
D(0) =
0 0.
’
9 i
0
-
&e
&J
,
(12)
. . . . . . . &d
. . . . . . . . . . . l/l,_
Y=Xl.
It will be shown later that all the eigenvalues of equations (11) are real. For convenience, the square root of 1: will be specified as follows : A=
1
I
if A: is positive
*
ifl
if 1: is negative.
It is easy to see that the solutions of equations (9) and (10) are P = XD(0) C, + XD-i(0) Cz,
w = AYD(B) C, - AYD-‘(0) Cz,
(13)
where C1 and C2 are two arbitrary constant column vectors. Now solutions’(l), (2), (3) and (13) must be joined together at the planes z = 0 or 0 = 0 and q = 0 or 0 = &. The appropriate boundary conditions are the continuity of pressure and velocity component normal to these planes. With the form of solutions given above, the imposition of these boundary conditions leads to at 8 = 0, P(r) + e,py) = XC1 + XC2,
EP”’ - ie, rc,p(jl)= AYC, - AYC2,
(14)
- EP”’ = AYDCl - AYD-’ C2,
(15)
at e = eE, P”) = XDCl + XD-’ C2, where D = D(e,) and E is a diagonal matrix zk, 0. .o . . . 0
97
SOUND TRANSMISSION IN DUCT BENDS
where K, = [CD* - a’(n - 1)”- Irzg*(k - 1)2]“2,
if o2 > 7r*(n- 1) + 7r2B2(k- 1)2
i[o’(n - l)* - 7r2/12(k- 1)2 - ~~]i/*,
if o* < 7r*(n- 1) + ~*/.?*(k- 1)2.
Equations (14) and (15) provide four vector equations for the four unknowns PC’), PC’), C1 and C2. By straightforward elimination of unknowns it is easy to find that the amplitudes of the reflected and transmitted waves are given by: P(‘) = Re,py ),
P(I) = Te,p’t) J’
(16)
where R and T are the reflection and transmission matrices. They can be computed according to the following formulas : R = R, + iR,, T = T(E + in, I), T = 2ZD-‘(EX R, = -Z[I
+ D-‘(EX + AY)-‘(EX
k, = Z[I - D-‘(EX + AY)-‘(EX
+ AY)-I,
- AY) D-l] X-l, - AY) D-‘1 Y-’ A-‘,
Z = [x-’ + Y-IA-’ E + D-‘(EX + AY)-‘(EX
- AY)D-'(X-l
-
Y-'A-’ E)]-‘.
(17)
(I is the identity matrix.)
In deriving expressions (17) effort has been made to isolate the indexj of the incident wave mode. Matrices T, R, and R, are independent ofj. Thus once they are computed they can be used for a wide range of j values. 3. SOME PROPERTIES OF THE GALERKIN
SOLUTION
The solution obtained by the Galerkin method in the previous section possesses certain properties which have been proved to be very valuable in the actual numerical computation of the transmission and reflection matrices. Moreover, some of the properties have also been found helpful in providing a physical and mathematical interpretation of the present solution. These, together with the relationship between the Galerkin solution and the classical solution by the method of separation of variables, will be examined below.
3.1.
THE EIGENVALUES 1; ARE REAL
From equations (11) the eigenvalues and eigenvectors are defined by -BAa, = ,ulai,
p(1= n:.
Since A and B are real and symmetric matrices, B-’ is also a real and symmetric matrix. Premultiplying the above equation by B-’ leads to -Aa, = n1 B-’ al.
(18)
The complex conjugate of equation (18) is -Aa: = &‘B-‘a:
(19)
(* denotes a complex conjugate). Premultiply equation (18) by a:’ and equation (19) by a:. (Superscript T indicates the transpose.) Upon subtracting one equation from the other it is easy to find that (c({- ,uy) a; B-la: Since a:BT’a: are real.
= 0.
# 0 it follows that pLI= pf so that the eigenvalues and eigenvectors of -BA
98
C. K. W. TAM
3.2. THE LIMITING CASE OF R --f m As R -+ a, the curved duct effectively becomes a straight duct. In this limit it is easy to A - -RE2. The eigenvalues of -BA show that A and B are diagonal matrices: B + RI, are therefore -R2 K: or Aj/R = +iKj. The ~c~(j= 1,2,. .) are, of course, the axial wave numbers of a straight duct. The Galerkin solution is of the form P N aje’J8-i0’. In the limit R + to, this becomes P N a, e * iK~S-iWf,where s = RB is the axial distance measured along the duct. Clearly the straight duct solution is recovered in this limit. Further, if nf is negative, the wave number is real and hence the solution represents a propagating mode. On the other hand if ,Jf is positive the wave number is imaginary. In this case the solution gives a non-propagating or cut-off mode. It is found that this association of the eigenvalues of -BA to propagating and cut-off modes is very useful. This point will be discussed further later. 3.3. RELATIONSHIP
BETWEEN
EIGENVALUES
Ai AND
SEPARATION
CONSTANTS
OF CLASSICAL
APPROACH
Previous investigations [l-S] on this problem invariably attempted to solve the equation of motion (the wave equation) by the method of separation of variables. According to this approach a solution in the form p mf(r)cos[n/?(k - l)y]eive-iwr is being sought. A closer examination of this solution indicates that it is closely related to the one obtained by the present method. The differences between the two solutions lie in that while the method of separation of variables expresses f(r) in terms of Bessel functions of order v the Galerkin method expands f(r) as a linear combination of independent basis functions: f(r) = x:m”_, C,COS[TT(WZ - l)(r - R)]. The separation constants Vi (sometimes called the angular wave numbers) are equal to the jLi, the eigenvalues of our matrix equation. The present method approximates the functionf(r) by using a finite number of terms thus leading to an approximate set of eigenvalues. However, it effectively avoids the eigenvalue equation involving complicated Bessel functions of non-integer order which has so far proved to be exceedingly difficult to deal with. 4. IMPLEMENTATION
OF THE GALERKIN
METHOD
AND RESULTS
In order to be able to compute the reflection and transmission matrices of equation (16) efficiently it is necessary to have a rather simple computation procedure which solves the eigenvalue problem (11). Fortunately, there is available a QR algorithm (see the book by Wilkinson [l 11) which can handle precisely this aspect of the computation. Standard Fortran subroutines for the determination of the eigenvalues and eigenvectors of a matrix can be obtained from several sources : e.g., EISPACK, the International Mathematical and Statistical Libraries. All the computations carried out in this work used the subroutines provided by IMSL. Aside from the eigenvalue problem the only other difficulty that one might encounter in the computation of R and T is in the operation of matrix inversion. This point has been carefully borne in mind when equation (17) was derived. Because some of the eigenvalues are large and positive, terms such as exp(& 6,) are exceedingly large. They are potential sources of problems in matrix inversion. As can be seen in equation (17) these terms are now systematically isolated and grouped into the matrix D-l. D-’ is a diagonal matrix with elements exp(-1,,0,) which can be computed readily without having to perform numerical matrix inversion. In this way possible matrix inversion problems are solved. AS has been shown the eigenvalues 1: are real. In the present computation they are ordered according to their magnitudes with 2: being the smallest. This is a natural way of ordering the eigenvalues since by increasing the size of the truncated matrices larger and larger eigenvalues are added to the sequence without disturbing the established order. Also physically the negative eigenvalues are associated with the propagating modes and they are the modes that are
03 02) 03 02) 03 02) 04 03) 04 03) 04 03)
-2.145409E ( 1G44967E -2.145204E (1 GM9668 -2.145199E (1 G44966JZ -
-
-7.673766E (4.193724E -7.669504E (4.1941878 -7.669096E (4.194224E
+6.327702E (6.8015558 +6.327264E (6.803787E +6.3272238 (6.803575E
-
-
-2.2733068 (-2.8273568 -2.272394E (-2.826639E -2.272392E (-2.826613E
-
(1.0616588 -
+2.078421E (l-0616248 +2.078984E (1.061656E +2.079041E
-1.2639668 (-2.585571E -1.263979E (-25848008 -1.263982E (-2.584720E
+7.944810E (--9.922793E +7.94468OE (-9.923255E +7.944650E (-9.923299E
-2.021378E (-1.9573538 -2.02 1449E (-1.957246E -2.021457E (-1.957234E -
04 04) 04 04) 04 04)
03 02) 03 02) 03 02)
02 03) 02 03) 02 03)
03 03) 03 03) 03 03)
02 02) 02 02) 02 02) -
+1.984058E (-2.249696E +1.983929E (-2.2497368 +I.9839158 (-2.249745E
-
-8.726403E (2.137794E -8.725031E (2.137893E -8.724921 E (2.1379OOE -
-1.238386E (-2.5663418 -1.2387OOE (-2.5663 14E -1.238734E (-2.566311E -
- 1.263966E (-2.585571E -1.263979E (-2.5848OOE -1.263982E (-2.584720E
-2.145409E ( 1+l4967E -2.145204E (1 Gl48OOE -2.145199E (1 M4966E -
03 03) 03 03) 03 03)
03 02) 03 02) 03 02)
03 02) 03 02) 03 02)
02 03) 02 03) 02 03)
03 02) 03 02) 03 02) -
-
+4*121327E (3.889103E i4121631E (3.8886048 +4.1216288 (3.888565~
-2.855467E (1.845766E -2.855227E (1.8460388 -2.8552 11E (1.8460578
-
-
-8.7264Q3E (2.137794E -8.72503 1E (2.137893E -8.724921 E (2.1379OOE -
+2.078421E (1.061624E +2.078984E (I.0616568 +2.079041E ( 1.061658E
-7.673766E (4.1937248 -7.669504E (4.1941878 -7.669096E (4.194224E
03 03) 03 03) 03 03)
02 02) 02 02) 02 02)
03 02) 03 02) 03 02)
03 02) 03 02) 03 02)
04 03) 04 03) 04 03)
-1.4370048 (6.401786E -1.436426E (6.4019188 -1.4364028 (6401924E
-
-
-
+1.984058E (-2.249696E +1.983929E (-2.249736~ +1.983915E (-2.2497458 +4.1213278 (3.8891038 +4.1216318 (3.8886048 +4.1216288 (3.888565E
-
-
-2.273306E (-2.8273568 -2.272394E (-2.826639E -2.2723928 (-2.826613~
+6.327702E (6.8015558 +6*327264E (6.803787E +6*327223E (6.803575E
02 02) 02 02) 02 02)
03 03) 03 03) 03 03)
03 03) 03 03) 03 03)
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
used. Only thefirst 5 x 5 submatrix
w = 15.0, 0, = a/2, b = 1.0, k = 2, R = 3.0. The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
04 05) 04 05) 04 05)
02 02) 02 02) 02 02)
-
-2.021378E (-1.957353E -2.021449E (-1.957246E -2.021457E (-1.957234E
02 02) 02 02) 02 02)
-
+1*303077E (-1.009861E +1.303071E (-lGO9866E +1.303070E (-1.009872E
is shown below)
TABLE 1 Transmi$sion matrix T (There are$ve propagating modesfor the parameters
04 04) 04 04) 04 04)
05 04) 05 04) 05 04)
+4.88298OE (2.8368378 +4,883418E (2.835938E +4.883455E (2~835895E -
-
-
-
-8.69345OE (1.881421E -8.738582E (1.881236E -8.741247E (1.88123lE
+6.565637E (-lW7883E +6.561418E (-1 G47079E +6.561144E (--1.047054E
+6.786596E (-6.105665E +6.781811E (-6.107183E +6*781503E (-6.107251E
-
+l.O15965E (1.7999OlE +l.O15623E (1.802748E +l.O156OlE (1.80291SE
-
-
+6444452E (-2.244243E +6.438010E (-2.2507078 +6.4375948 (-2.25114OE
+8.15968OE (-1.247656E +8.1585108 (-1.250408E +8-l 58445E (-1.2505728
-5WO768E (-1.508793E -5WO774E (-1.512977E -5+300774E (-1.513187E
+4.769656E (2.770999E i4770084E (2,770122E +4,770120E (2.77008OE -
03 04) 03 04) 03 04)
04 05) 04 05) 04 05)
04 04) 04 04) 04 04)
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
-
-5443959E (-9822010E -5449999E (-9.821526E -5.450338E (-9.821465E
-
+6443038E (-9.300689E +6441785E (-9.26903 1E +6441669E (-9.267227E -
-4996537E (2.235231E -4.9965398 (2.235711E -4.996539E (2.235754E
+7.5482578 (-1.1541668 +7.547 174E (-1.156712E +7,547114E (-1.156864E -
-7.855391E (1.70005OE -7.896172E (1.699883E -7.898581E (1.699878E -
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
05 04) 05 04) 05 04)
-
-
-5.004256E (-7.1759OOE -5.004258E (-7.172510E -5.004258E (-7.172224E +8.900146E (-1.259086E +8.900764E (-1.258989E +8.90077OE (-1.258977E
-
-
+5.055236E (-1.760457E +5.050183E (-1.765527E +5+M9857E (-1.765867E +5.463521E (-7.8867328 +5*462459E (-7.859887E +5.462361E (-7.858357E
-
04 03) 04 03) 04 03)
01 04) 01 04) 01 04)
04 05) 04 05) 04 05)
04 05) 04 05) 04 05)
04 05) 04 05) 04 05)
used)
+5.03077OE (-8.029163E +5.027537E (--8.023005E +5.0273278 (-8.022809E
-5.001183E (-3,088994E -5.001183E (-3.089944E -5+M1183E (-3.08995OE
+5.990322E (-8.47439OE +5.990737E (-8.473738E +5.990742E (-8.473659E
-3.107061E (-5605771E -3.110508E (-5605495E -3.110702E (-5.60546OE
$5.363981 E (9.5029228 +5.362178E (9.517953E +5.362061E (9.518833E
+3,499956E (-3.148789E +3.497488E (-3.149572E f3.497330E (-3.1496078
-
-
-
-
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
___ - 04 - 04) - 04 - 04) - 04 - 04)
The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
04 04) 04 04) 04 04)
04 04 04 04) 04 04)
01 05) 01 05) 01 05)
-
-5~000759E (-8.849512E -5@00763E (-8.874411E -5.000763E (-8.875728E
TABLE 2 Matrix RI (see Table I .for values qf various parameters
06 05) 06 05) 06 05) 05 05) 05 05) 05 05)
-7.1443358 (-4.476370E -7.1388558 (-4.473493E -7.13868lE (-4.473307E -4.162767E (-4.627017E -4.163802E (-4.6237548 -4.163848E (46235458
05 02) 05 02) 05 02) 06 05) 06 05) 06 05)
-
-1.053117E (-3.489396E -1.0560378 (-3.489392E -1.056184E (-3.4893928 -8.7084678 (-5.695346E -8.727677E (-5.694529E -8.728827E (-5.694484E -1.566451E (-4.49814OE -1.570963E (-4.493644E -1.571265E (-4.493354E +1.2563078 (-7.0912978 +1.258294E (-7.088915E +1.258410E (-7.088759E
05 05) 05 05) 05 05)
06 05) 06 05) 06 05)
05 05) 05 05) 05 05)
+1*93412OE (-3.329156E +1.9335088 (-3.329454E +1.9334788 (-3.329479E -
06 05) 06 05) 06 05)
-
-7.0176OOE (-4.861431E -6.993713E (4.860486E -6.992352E (-4.860399E -7.4109498 (4.107602E -7.410584E (4.1121598 -7.4105388 (4.112415E
05 05) 05 05) 05 05)
05 06) 05 06) 05 06) 06 05) 06 05) 06 05) 05 02) 05 02) 05 02)
+1.282730E (5.927087E +1.282603E (5.957857E +1.2826OOE (5.9596748 -8.708468E (-5.6953468 -8.7276778 (-5.694529E -8.7288278 (-5.694484E +1.686537E (-3.7752378 +1.6868998 (-3.775235E +1.686932E (-3.775235E 05 02) 05 02) 05 02) 04 05) 04 05) 04 05)
-6.385103E (4445204E -6.382086E (4445203E -6.381832E (-4445203E -1.120332E (-7.919334E -1.120246E (-7.919883E -1.120236E (-7.919889E
-
06 05) 06 05) 06 05)
-
-7.0176OOE (4861431E -6.993713E (-4860486E -6.992352E (4.860399E
06 05) 06 05) 06 05) 06 05) 06 05) 06 05)
-
-7.1443358 (447637OE -7.138856E (4.4734938 -7.138681E (4473307E -1.566451E (449814oE -1.570963E (-4.493644E -1.571265E (-4.4933548
-
04 05) 04 05) 04 05) 4.083716E - 05 (-6.608544E - 02) -4.084972E - 05 (-6608543E - 02) 4084980E - 05 (-6608543E - 02)
-1.120332E (-7.919334E -1.1202468 (-7.919883E -1.120236E (-7.919889E
-4.162767E - 05 (4.627017E - 05) -4.163802E - 05 (4.6237558 - 05) -4.163848E - 05 (4.6235458 - 05) +1*256307E - 05 (-7.0912978 - 05) +1.258294E - 05 (-7.088915E - 05) +1.258410E - 05 (-7.0887598 - 05) -7*41095OE - 05 (4.107602E - 05) -7.410584E - 05 (4.112159E - 05) -7.410538E - 05 (4.112415E - 05)
The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
-
06 02) 06 02) 06 02) 05 05) 05 05) 05 05) 05 06) 05 06) 05 06)
-6.0334888 (-3408421E -6.050463E (-3408418E -6.0513618 (-3408418E +1.934120E (-3.329156E +1.933508E (-3.329454E +I .933478E (-3.329479E +I .28273OE (5.927087E +1.2826048 (5*957857E +1+2826OOE(5.959675E -
TABLE 3
Matrix R, (see Table 1 for values of various parameters used)
102
C. K. W. TAM TABLE 4
Eigenvalues (A:) w = 40x, /I = 1.O, k = 11, R = 1 (only thefirst 40 eigenvalues are tabulated below) N= 100
N=50
__1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
-6.003698E -5.280253E -4.8755768 -4.572737E -4.3 10502E -4.079362E -3.869206E -3.676555E -3.497749E -3.330858E -3.174172E -3.026407E -2.886661E -2.753989E -2.62788OE -2.507594E -2.392832E -2.282982E -2.177865E -2.0769538 -3.980124E -1.886885E -1.797106E -1.710183E -1.625543E -1.541626E -1.459617E -1.406372E -1.34694OE -1.25248OE -1.148122E -1.0371OOE -9.202061E -7,978 184E -6.701889E -5.374953E -3.99876OE -2.574394E -1.102737E +4 154960E
+ 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 -t 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 03 + 03 + 03 + 03 -t- 03 -t 03 + 03 + 02
-6.003699E -5.280254E -4.875576E -4.572737E -4.3 10502E -4.079362E -3.869206E -3.676555E -3.497749E -3.330858E -3.174172E -3.026407E -2.886661E -2.753989E -2.627880E -2.507594E -2.39283 1E -2.282982E -2.177865E -2.076953E -1.980124E -1.886885E -1.797106E -1.710183E -1.625543E -1.541626E -1.459616E -14063728 -1.346941E -1.25248OE -1.148122E -1.037lOOE -9.20206OE -7.978 I80E -6.791885E -5.374946E -3.998753E -2.574386E -1.102736E +4.154727E
+ 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 -I- 04 -I- 04 + 04 f 04 + 04 + 04 + 04 + 04 + 04 + 04 + 03 + 03 + 03 + 03 + 03 + 03 + 03 + 02
of practical interest. With the above ordering of eigenvalues the relevant information about these modes are contained in the upper left-hand submatrices of R and T. Thus only these submatrices need be considered in actual computation. Tables 1, 2 and 3 represent some typical results of a routine computation. The numerical values of the parameters used are w = 15.0, p = 1.0, k = 2, f3, = 71/2 and R = 3.0. In this case there are five negative eigenvalues and hence five propagating modes. The first 5 x 5 submatrices of the matrices T, RI and Rz are given in the tables. Three sets of values with sizes of the matrices truncated to N = 10, 20 and 25 are exhibited. It can be seen that there is a rapid numerical convergence using the present computation procedure. Table 4 provides an example of numerical results for computations involving high frequency and high mode number incident wave. The parameters are w = 4071, ,9 = 1.O, k = 11 and R = 1
SOUND TRANSMISSION IN DUCT BENDS
103
25 -
w
Figure 3. Eigenvalues of propagating
modes. k = 2; fi= 1.0; R = 3.0.
(very sharp bend). There are thirty-nine negative eigenvalues and hence thirty-nine propagating modes (with the same k value). Given in the first column of this table are the first forty eigenvalues computed by using truncated matrices of 50 x 50. They are to be compared with the first forty eigenvalues obtained by using matrices of 100 x 100 shown in the second column. It is clear from such a comparison that the negative eigenvalues are accurate to at least six significant figures. This accuracy is sufficient for most applications. Figure 3 provides an example illustrating the dependence of Im(&), or the angular wave number, on the frequency, o. Only the eigenvalues of the propagating modes are shown. As frequency increases there are more and more propagating modes just as in the case of a straight duct. It is to be noted that even at moderately large values of o the angular wave numbers of the propagating modes are quite large. This implies (verified by actual computation) that the phases ofthe transmission and reflection coefficients will oscillate very rapidly between positive and negative values even with small changes in 19~or CD.Thus in considering the transmission of sound through a duct system good accuracy must be retained in joining the solutions together at the two ends of a curved bend so as to ensure correct matching of phases. 5. DISCUSSION
To make the above computation procedure practical it is necessary to truncate the matrices to a finite size of N x N. Naturally, the larger the N is the more accurate the computed solution will be. For some problems, however, only a reasonably accurate solution is all that is required. In these cases it is meaningful to ask what is the smallest N for which such a solution could be expected. To estimate the value of Nminit is worthwhile to note that the eigenvalues of nt are real. Further the negative eigenvalues give rise to propagating modes while the positive eigenvalues provide the cutoff modes. For a bend with a not too small turning angle dE, the cut-off modes play only a very insignificant role in the transmission of sound. They are rapidly attenuated spatially. This observation indicates, therefore, that Nmin is an integer somewhat larger than the total number of negative eigenvalues. This is confirmed computationally. An illustration of this is given by the very accurate results of the first column of Table 4 where N = 50 (39 negative eigenvalues) is used.
104
C. K. W.
TAM
One interesting
phenomenon observed in our numerical experiments testing the present procedure is that the transmission matrix T is symmetric (see Table 1). From the second of equations (17) it follows that T is also symmetric and that the index j of the incident wave appears only in the diagonal elements of T. The physical implication is that a reciprocal relationship between incident and transmitted wave modes exists. That is to say that the transmission coefficient for the mth modes due to an incident wave of the nth mode is the same as the transmission coefficient of the nth mode due to an incident wave of the mth mode. This
interesting property does not seem to have been noticed before. In the work reported here the computation procedure has been developed primarily for rectangular ducts with rigid walls. However, it is straightforward to see that with only simple modifications the same procedure could be used for ducts with soft walls. Moreover, by using a toroidal co-ordinate system instead of a cylindrical co-ordinate system the present method is applicable to the study of sound transmission in circular bends of circular ducts. In this case the wave equation is non-separable so that the classical method of separation of variables cannot be used. It is envisaged that computation results involving soft walled rectangular ducts and circular ducts will be reported in a subsequent paper. ACKNOWLEDGMENT This work was supported in part by the National Science Foundation under Grant GK35790 and the National Aeronautical and Space Administration Langley Research Center Grant NSG 1021. The author wishes to thank Barbara Schreur for her assistance in computer programming. REFERENCES 1. P. E. KRASNUSHKIN 1945 Uch. Zap. Mosk. Gos. Univ. No. 75, Bk. 2, Pt. 2, 9-27.
On waves in curved tubes. 2. A. G. SEVESHNIKOV1958 Radiotekhnika i Eiektrouika 3, 641-648. Waves in bent tubes. 3. F. E. GRIGOR’YAN 1969 Akustische Zeitschrift 14, 376-384. (English translation: Soviet Physics Acoustics 14, 315-321, 1969). Theory of sound wave propagation in curvilinear wave guides. 4. W. C. OSBORNE1968 National College for Heating, Ventilating, Refrigeration and Fan Engineering Technical Memo. No. 5. The propagation of sound waves round a bend of rectangular section. 5. W. ROSTAFINSKI1972 Journalof the AcousticalSociety ofAmerica52,141 I-1420. On propagation of long waves in curved ducts. 6. W. ROSTAFINSKI1974 Journal of the Acoustical Society of America 56,l l-15. Analysis of propagation waves of acoustic frequencies in curved ducts. 7. A. CUMMINGS1974 Journal of Sound and Vibration 35, 451-477. Sound transmission in curved duct bends. 8. M. K. MYERS and P. MUNGUR 1975 American Institute of Aeronautics and Astronautics Paper 75-497. Sound propagation in curved ducts. 9. B. A. FINLAY~~N1972 The Method of Weighted Residuals and Variational Principles. New York: Academic Press. 10. W. EVERSMAN,E. L. COOK and R. J. BECKEMEYER1975 Journal of Sound and Vibration 38, 105-123. A method of weighted residuals for the investigation of sound transmission in nonuniform ducts without how. 11. J. H. WILKINSON 1965 The Algebraic Eigenvalue Problem. The Oxford University Press.
91-104
A STUDY OF SOUND TRANSMISSION DUCT BENDS BY THE GALERKIN
IN CURVED METHOD
C. K. W. TAM Department of Mathematics, Florida State University, Tallahassee, Florida 32306, U.S.A.
(Received 19 June 1975, and in revisedform
8 August 1975)
A computational procedure based on the Galerkin method is developed to study the sound transmission and reflection characteristics of circularly curved bends of rectangular ducts. This procedure is adaptable to routine computer calculation. An important component of the computer program consists of a QR iterative algorithm which solves the matrix eigenvalue problem generated by the Galerkin method. The present procedure produced very satisfactory convergent results when applied even to cases where the incident wave has high frequency and high wave mode numbers. Some properties of the Galerkin solution are investigated. Its relationship to the classical solution by the method of separation of variables is discussed. 1. INTRODUCTION
This paper is concerned with the computational aspect of the transmission of sound through circularly curved bends of a duct system with rectangular cross-section. Curved bends are unavoidable features of most practical duct systems. An efficient procedure to compute their acoustical characteristics should, therefore, be valuable to engineers and architects. The propagation of sound waves in curved ducts is characterized by wave patterns which differ in many important respects from those in straight ducts. In a duct system consisting of straight ducts linked together by curved bends the continuous reflection of sound waves from the curved boundary of a bend causes a partial reflection of any incident wave from the adjoining straight ducts. This reduces the sound transmission capability of the duct system as compared with that of a single straight duct. The physics of the problem is rather simple and can easily be understood. The main problem from a practical standpoint is computational difficulty. For a circular bend of a rectangular duct it was recognized by Krasnushkin [l] and others, some years ago, that the equations of motion were separable with respect to a suitably chosen cylindrical co-ordinate system. However, for a bend of finite length this solution must be joined to the incident, transmitted and reflected wave solutions of the adjoining straight ducts. To avoid these problems Krasnushkin and later Seveshnikov [2], Grigor’yan [3] and others restricted themselves to the consideration of infinitely long bends. This is, of course, an idealization and it corresponds approximately to that of a helical wave guide. Physically, their findings could therefore provide only qualitative information for application to practical duct systems. Seveshnikov [2] generalized the problem considered by Krasnushkin [I]. He investigated theoretically acoustic wave propagation in arbitrarily deformed wave guides. On starting from the acoustic wave equation written for a general polar co-ordinate system he studied solutions which would correspond to a helical elbow bend of a circular section pipe. Unfortunately, the validity of his results is restricted to ducts with bends of very large radius of 91
92
C. K. W. TAM
curvature to diameter ratios. Osborne [4] investigated theoretically the separable solution of sound wave propagation in rectangular section curved ducts. The radial equation leads to Bessel functions of non-integer order which are exceedingly difficult to compute. He limited his studies to cases of very large inner and outer bend radii and very small inner radius. Recently works on the subject of sound transmission in curved rectangular ducts were published by Rostafinski [5, 61, Cummings [7] and Myers and Mungur [8]. Their approach follows very much the method of separation of variables used by previous investigators. However, Rostafinski and Cummings did make an effort to join the solution in the curved part of the duct to the propagating wave solutions in the straight portion of the duct system. Again the major problem of having to deal with Bessel functions of non-integer order forced Rostafinski to consider only very low frequency sound and Cummings to study only the propagation of the lowest (plane) wave mode. In his more recent publication Rostafinski [6] tried to avoid the problem of computing Bessel functions of non-integer order by treating the problem in an inverse manner. In this way he was able to produce some numerical results including velocity distributions across the duct. Unfortunately, this method does not guarantee that the full set of eigenfunctions are being used. Thus, it is possible that the computed results are incomplete. Moreover, this method, even if it is free of the above difficulty, is very cumbersome to use in a practical situation because it does not deal with the given physical problem directly. The primary objective of this investigation is to develop a direct, reliable computational procedure which can be used to determine the sound transmission and reflection characteristics of circular bends of rectangular ducts. The procedure is to be adaptable to routine computer calculations so that it can provide useful information for design or be used as a guide for other engineering purposes where sound transmission is a factor to be taken into consideration. For this purpose, the method of separation of variables which has been used by all previous investigators is found to be inadequate. In this paper the method of Galerkin modified to suit the present problem is adopted instead. The Galerkin method has been used by many authors to provide solutions to various engineering problems. A description of this method can be found in books (e.g., that of Finlayson [9]) and will not be repeated here. In duct acoustics, this method has recently been employed successfully by Eversman, Cook and Beckemeyer [lo] to study the propagation of sound in ducts with variable cross-section. For the problem under consideration, the Galerkin method, as will be shown in section 2, reduces the wave equation to a large system of linear ordinary differential equations with constant coefficients. It is well known that the solution of these equations can be constructed in terms of the eigenvectors and eigenvalues of the coefficient matrix. Thus, the success of the Galerkin method hinges on how efficiently the matrix eigenvalues problem can be solved numerically. Fortunately, there is available a stable iterative QR algorithm (see the book by Wilkinson [I 11) which is particularly effective in handling this aspect of the computation. The combined Galerkin method-QR algorithm has now been programmed into a special subroutine at the Florida State University Computer Center. Through direct numerical experiments using the Center’s CDC 6500 Computer, it is found that the procedure is well suited for the goals stated. The cases tested included high frequency, high order wave modes incident on a sharp circular bend. These tests are believed to cover not only most of the practical cases but also a wide range of parameters for which the method of separation of variables used previously failed to produce fruitful numerical results. Some properties of the solution obtained by the Galerkin method will be presented in section 3. These properties are found to be very helpful in interpreting the solution physically and mathematically. Further, they provide certain information which has proved to be valuable in actual numerical computations. Some experience gained in applying the present computation procedure will be discussed at the end of this paper.
SOUND TRANSMISSION IN DUCT BENDS
93
2. FORMULATION
Consider the transmission of sound through a circularly curved section of a duct system as shown in Figure 1. Here it is assumed that the cross-section of the duct is the same on both
I
Figure 1. Curved duct geometry.
sides of the bend. Without loss of generality the problem is conceived to consist of an incident wave with mode numbers (j,k) propagating in the positive z-direction along region (1) of the duct system. On reaching the curved portion of the duct, namely region (2), part of the wave would be reflected. The reflected wave propagates in the negative z-direction away from the x-y plane. The other part of the incident wave is transmitted into region (3) and radiates away from the c-u plane. Let p, U, u and w denote the pressure and the velocity components in the x, y and z, or 5, y and q (see Figure 2), directions, respectively. In terms of
7
Figure 2. Co-ordinate systems used.
dimensionless variables with length scale b (width of duct in the x-direction), velocity scale c (the speed of sound) time scale b/c and pressure scale pc2 (p is the gas density), the incident, reflected and transmitted waves in regions (1) and (3) can be written as follows : the incident wave (j, k)t mode : P (i)=p~)cos[lc(j2.P = i{Tc(jvti) = i{$(k
I)x]cos[nj?(k-
l)/~}p:~)sin[7c(j- l)/o}pj”
w(*) = (~~/o)p~~)cos[dj-
cos[+-
l)y]e’(KJL-W’), 1) x] cos[$(k l)x]sin[$(k
- l)y] e’(KJ’-“f), - l)y]e’(KJz-Or,r),
1)x] cos[nfi(k - I)y]ei(“Jz-‘Ot),
(1)
t Throughout this paper the wave modes are numbered in such a way that the plane wave mode is represented byj=l,k=l.
94
C. K. W. TAM
where ~~ = [co2 - n'(j- 1)”- rc2(k- l)*/?*]i/*, p,w is the amplitude of the incident wave and B = b/h is the aspect ratio of the duct; the reflected waves: P (r) =f&p;’
cos[7r(m -
1) X]cos[rcfi(k - l)y] e-i(Kmr+~t),
u(‘) =mzl i{rc(m - l)/o}pg) sin[rr(m - 1) X] cos[r$(k - l)y] e-i(Kmz+wt), n(r) = $i i{@(k - 1)/o] P,(‘) cos[lr(m
-
1) x] sin[ nj3(k - 1) y] e--i(Kmz+of),
WC’=) $,-(Km/cO)P(m’)cos[rc(rr?- 1) x] cos[rcj?(k - l)y] e-i(K~z+of),
(2)
where if w2 > n’(m - l)* + x2 /I*& - l)z,
K, = [W’- 7?()72- 1)2 - 7r2jY2(k- 1)2]1/2; ( i[n*(m - l)* + 7r2j12(k - 1)2 - ~0~1~‘~;
if o2 < 7r*(m- 1)2 + 7c2j12(k - 1)2;
the transmitted waves: P (f) =~lpg)
cos[7c (m - 1) [] cos[@(k - l)y] e’(xmn--of),
u(‘) =zii(a(m
- 1)/w} py’ sin[7r(m - 1) c] cos[n/?(k - 1) y] ei(emq-a)r),
- 1) 51sin[ a/?+ - 1) y] ei(Km*-aof), u(‘) =mzl i{rr(k - 1) /3/o}ps’ cos[ ~~(112
w(‘)=~l(K”/4P:)
cos[ 7t(m - 1) [] cos[ aj?(k - 1) y] ei(~mU-or).
(3)
Equations (I), (2) and (3) are, of course, solutions of the equations of motion,
ap _67+ r*v=o,
$&VP.
(4)
They satisfy the rigid wall boundary conditions at the surface of the duct. In addition, the radiation condition has been invoked in selecting the form of the reflected and transmitted waves. The amplitudes (complex in general) p$) and pg’ are unknown at this stage. The determination of their relationship to pj*) is the central problem of this work. In the curved part of the duct system i.e., region (2) of Figure 1 and Figure 2, the equations of motion in cylindrical co-ordinates are
ap au i aw u -+-+--+-+-~~, at ar r ae r ;+g=o,
au ay
au ap -+--0, at ay
fY+!&o.
(5)
A solution of equations (5) will now be constructed by the Gale&in method. According to this method the dependent variables are to be represented by a complete set of basis functions with arbitrary amplitudes. Here, the following basis functions will be used : p = J,
p,(O) cos[w(m - 1) (r - R)] cos[nfi (k - 1) y] e-‘“‘,
u = mz, i{lr(m - l)/o}p,(8)
sin[rc(m - l)(r - R)] cos[r$(k - 1) y] e-i”,
SOUND TRANSMISSION
v
=mzl i(7c(k-
1) /?/o}pJe)
95
IN DUCT BENDS
cos[rr(m - l)(r - R)] sin[7$(k - l)y] e-‘“‘,
w = $, (i/w) w,(8) cos[rr(m - 1) (r - R)] cos[r$(k - 1) y] e-rat.
(6)
Here R Q r < R + I, 0 < y < l/p, R being the inner radius of curvature of the curved duct. In principle, there is very little restriction as to what type of function could be used as basis functions so long as the functions form a complete set preferably satisfying certain boundary conditions. However, it is well known that for a given problem some basis functions work better than others. For this reason, it is always advantageous to use basis functions that are physically or mathematically natural to the problem. The basis functions of equations (6) are chosen for two reasons. First, they are simple and are closely related to the duct modes of a straight duct. Secondly, at f3= 0 and 0 = eE these functions match automatically with the incident, reflected and transmitted wave modes. This latter property will facilitate the process of joining these solutions together at these two planes. The two sets of functions Pm(B)and w,(e) are as yet arbitrary. They are to be chosen such that the continuity and &momentum equations of equations (5) are satisfied approximately. This is done by first substituting equations (6) into the first and fourth of equations (5) to obtain the residuals. Then the residuals are forced to become zero by requiring them to be orthogonal to a complete set of functions. For this purpose the following complete set of functions will be used : (7) r cos[n(n - 1) (r - R)], n=l,2,3,... The residuals are multiplied by each of the functions (7) and set equal to zero after integration over r from R to R + 1. In this way two infinite sets of equations for pm(e) and w,,,(e)are found. In practice these systems of equations and also the summations in equations (2), (3) and (6) are all truncated at a large integer value : say, M = n = N. The resulting finite system of equations and unknowns are hence rendered solvable. On carrying out the above steps it is easy to find that&@ and w,(e) are given by the solution of the following systems of ordinary differential equations with constant coefficients. dw,/de = 5 A,,P,,
dp,W’ = - 5 B,,, w,,
(n
=
1,2,3,...N)
where ;
&m =
[(-l).+,-
l]
[(m - 1)2 + (n - 1)2] [(m - I)2 - (n - l)‘]’ ’
n#m
_
I
n=rn
;+R, [o2
A,,,, =
(8)
m=l
m=l
_
x2
B2(k
_
1)21
B”,
_4(m- II2(n- II2rel)m+n- 11) [(m - 1)2 - (n - 1)2]2
[02-~2j?2(k-1)2-~2(m-l)2](1/2+R),
n#m n=m.
To facilitate the solution of equations (8) the following matrix notation will be introduced A = matrix [A,,], B = matrix [B,,,], P = column vector with elements p#), w = column vector with elements w,(e), p(‘) = column vector with elements p’,), p@)= column vector with elements pz), e, = column vector whosejth element is unity and all the other elements are zero.
96
c. K. w. TAM
It is important to note that the matrices A and B are real and symmetric. In matrix notation equations (8) can be rewritten as dP/dQ = -Bw.
dw/d0 = AP,
(9)
On eliminating w from equations (9), a single equation for P is found: d2 P/de2 = -BAP.
(10)
A general solution of equation (10) can easily be constructed in terms of the eigenvalues and eigenvectors of the matrix -BA. Suppose A:, A:. . . A: are the eigenvalues and CQ,a,, . . . aN are the corresponding eigenvectors of -BA: i.e., -BAa, = 1: a,.
(11)
The following matrices then can be formed : X = [al, a2,. . . aN], a matrix whose ith column is ai, l/A1
0
0
0
l/l2
0
0
l/A,
A=0
. . 0
r
eAle 0 0. . . . 0
-
0
:
D(0) =
0 0.
’
9 i
0
-
&e
&J
,
(12)
. . . . . . . &d
. . . . . . . . . . . l/l,_
Y=Xl.
It will be shown later that all the eigenvalues of equations (11) are real. For convenience, the square root of 1: will be specified as follows : A=
1
I
if A: is positive
*
ifl
if 1: is negative.
It is easy to see that the solutions of equations (9) and (10) are P = XD(0) C, + XD-i(0) Cz,
w = AYD(B) C, - AYD-‘(0) Cz,
(13)
where C1 and C2 are two arbitrary constant column vectors. Now solutions’(l), (2), (3) and (13) must be joined together at the planes z = 0 or 0 = 0 and q = 0 or 0 = &. The appropriate boundary conditions are the continuity of pressure and velocity component normal to these planes. With the form of solutions given above, the imposition of these boundary conditions leads to at 8 = 0, P(r) + e,py) = XC1 + XC2,
EP”’ - ie, rc,p(jl)= AYC, - AYC2,
(14)
- EP”’ = AYDCl - AYD-’ C2,
(15)
at e = eE, P”) = XDCl + XD-’ C2, where D = D(e,) and E is a diagonal matrix zk, 0. .o . . . 0
97
SOUND TRANSMISSION IN DUCT BENDS
where K, = [CD* - a’(n - 1)”- Irzg*(k - 1)2]“2,
if o2 > 7r*(n- 1) + 7r2B2(k- 1)2
i[o’(n - l)* - 7r2/12(k- 1)2 - ~~]i/*,
if o* < 7r*(n- 1) + ~*/.?*(k- 1)2.
Equations (14) and (15) provide four vector equations for the four unknowns PC’), PC’), C1 and C2. By straightforward elimination of unknowns it is easy to find that the amplitudes of the reflected and transmitted waves are given by: P(‘) = Re,py ),
P(I) = Te,p’t) J’
(16)
where R and T are the reflection and transmission matrices. They can be computed according to the following formulas : R = R, + iR,, T = T(E + in, I), T = 2ZD-‘(EX R, = -Z[I
+ D-‘(EX + AY)-‘(EX
k, = Z[I - D-‘(EX + AY)-‘(EX
+ AY)-I,
- AY) D-l] X-l, - AY) D-‘1 Y-’ A-‘,
Z = [x-’ + Y-IA-’ E + D-‘(EX + AY)-‘(EX
- AY)D-'(X-l
-
Y-'A-’ E)]-‘.
(17)
(I is the identity matrix.)
In deriving expressions (17) effort has been made to isolate the indexj of the incident wave mode. Matrices T, R, and R, are independent ofj. Thus once they are computed they can be used for a wide range of j values. 3. SOME PROPERTIES OF THE GALERKIN
SOLUTION
The solution obtained by the Galerkin method in the previous section possesses certain properties which have been proved to be very valuable in the actual numerical computation of the transmission and reflection matrices. Moreover, some of the properties have also been found helpful in providing a physical and mathematical interpretation of the present solution. These, together with the relationship between the Galerkin solution and the classical solution by the method of separation of variables, will be examined below.
3.1.
THE EIGENVALUES 1; ARE REAL
From equations (11) the eigenvalues and eigenvectors are defined by -BAa, = ,ulai,
p(1= n:.
Since A and B are real and symmetric matrices, B-’ is also a real and symmetric matrix. Premultiplying the above equation by B-’ leads to -Aa, = n1 B-’ al.
(18)
The complex conjugate of equation (18) is -Aa: = &‘B-‘a:
(19)
(* denotes a complex conjugate). Premultiply equation (18) by a:’ and equation (19) by a:. (Superscript T indicates the transpose.) Upon subtracting one equation from the other it is easy to find that (c({- ,uy) a; B-la: Since a:BT’a: are real.
= 0.
# 0 it follows that pLI= pf so that the eigenvalues and eigenvectors of -BA
98
C. K. W. TAM
3.2. THE LIMITING CASE OF R --f m As R -+ a, the curved duct effectively becomes a straight duct. In this limit it is easy to A - -RE2. The eigenvalues of -BA show that A and B are diagonal matrices: B + RI, are therefore -R2 K: or Aj/R = +iKj. The ~c~(j= 1,2,. .) are, of course, the axial wave numbers of a straight duct. The Galerkin solution is of the form P N aje’J8-i0’. In the limit R + to, this becomes P N a, e * iK~S-iWf,where s = RB is the axial distance measured along the duct. Clearly the straight duct solution is recovered in this limit. Further, if nf is negative, the wave number is real and hence the solution represents a propagating mode. On the other hand if ,Jf is positive the wave number is imaginary. In this case the solution gives a non-propagating or cut-off mode. It is found that this association of the eigenvalues of -BA to propagating and cut-off modes is very useful. This point will be discussed further later. 3.3. RELATIONSHIP
BETWEEN
EIGENVALUES
Ai AND
SEPARATION
CONSTANTS
OF CLASSICAL
APPROACH
Previous investigations [l-S] on this problem invariably attempted to solve the equation of motion (the wave equation) by the method of separation of variables. According to this approach a solution in the form p mf(r)cos[n/?(k - l)y]eive-iwr is being sought. A closer examination of this solution indicates that it is closely related to the one obtained by the present method. The differences between the two solutions lie in that while the method of separation of variables expresses f(r) in terms of Bessel functions of order v the Galerkin method expands f(r) as a linear combination of independent basis functions: f(r) = x:m”_, C,COS[TT(WZ - l)(r - R)]. The separation constants Vi (sometimes called the angular wave numbers) are equal to the jLi, the eigenvalues of our matrix equation. The present method approximates the functionf(r) by using a finite number of terms thus leading to an approximate set of eigenvalues. However, it effectively avoids the eigenvalue equation involving complicated Bessel functions of non-integer order which has so far proved to be exceedingly difficult to deal with. 4. IMPLEMENTATION
OF THE GALERKIN
METHOD
AND RESULTS
In order to be able to compute the reflection and transmission matrices of equation (16) efficiently it is necessary to have a rather simple computation procedure which solves the eigenvalue problem (11). Fortunately, there is available a QR algorithm (see the book by Wilkinson [l 11) which can handle precisely this aspect of the computation. Standard Fortran subroutines for the determination of the eigenvalues and eigenvectors of a matrix can be obtained from several sources : e.g., EISPACK, the International Mathematical and Statistical Libraries. All the computations carried out in this work used the subroutines provided by IMSL. Aside from the eigenvalue problem the only other difficulty that one might encounter in the computation of R and T is in the operation of matrix inversion. This point has been carefully borne in mind when equation (17) was derived. Because some of the eigenvalues are large and positive, terms such as exp(& 6,) are exceedingly large. They are potential sources of problems in matrix inversion. As can be seen in equation (17) these terms are now systematically isolated and grouped into the matrix D-l. D-’ is a diagonal matrix with elements exp(-1,,0,) which can be computed readily without having to perform numerical matrix inversion. In this way possible matrix inversion problems are solved. AS has been shown the eigenvalues 1: are real. In the present computation they are ordered according to their magnitudes with 2: being the smallest. This is a natural way of ordering the eigenvalues since by increasing the size of the truncated matrices larger and larger eigenvalues are added to the sequence without disturbing the established order. Also physically the negative eigenvalues are associated with the propagating modes and they are the modes that are
03 02) 03 02) 03 02) 04 03) 04 03) 04 03)
-2.145409E ( 1G44967E -2.145204E (1 GM9668 -2.145199E (1 G44966JZ -
-
-7.673766E (4.193724E -7.669504E (4.1941878 -7.669096E (4.194224E
+6.327702E (6.8015558 +6.327264E (6.803787E +6.3272238 (6.803575E
-
-
-2.2733068 (-2.8273568 -2.272394E (-2.826639E -2.272392E (-2.826613E
-
(1.0616588 -
+2.078421E (l-0616248 +2.078984E (1.061656E +2.079041E
-1.2639668 (-2.585571E -1.263979E (-25848008 -1.263982E (-2.584720E
+7.944810E (--9.922793E +7.94468OE (-9.923255E +7.944650E (-9.923299E
-2.021378E (-1.9573538 -2.02 1449E (-1.957246E -2.021457E (-1.957234E -
04 04) 04 04) 04 04)
03 02) 03 02) 03 02)
02 03) 02 03) 02 03)
03 03) 03 03) 03 03)
02 02) 02 02) 02 02) -
+1.984058E (-2.249696E +1.983929E (-2.2497368 +I.9839158 (-2.249745E
-
-8.726403E (2.137794E -8.725031E (2.137893E -8.724921 E (2.1379OOE -
-1.238386E (-2.5663418 -1.2387OOE (-2.5663 14E -1.238734E (-2.566311E -
- 1.263966E (-2.585571E -1.263979E (-2.5848OOE -1.263982E (-2.584720E
-2.145409E ( 1+l4967E -2.145204E (1 Gl48OOE -2.145199E (1 M4966E -
03 03) 03 03) 03 03)
03 02) 03 02) 03 02)
03 02) 03 02) 03 02)
02 03) 02 03) 02 03)
03 02) 03 02) 03 02) -
-
+4*121327E (3.889103E i4121631E (3.8886048 +4.1216288 (3.888565~
-2.855467E (1.845766E -2.855227E (1.8460388 -2.8552 11E (1.8460578
-
-
-8.7264Q3E (2.137794E -8.72503 1E (2.137893E -8.724921 E (2.1379OOE -
+2.078421E (1.061624E +2.078984E (I.0616568 +2.079041E ( 1.061658E
-7.673766E (4.1937248 -7.669504E (4.1941878 -7.669096E (4.194224E
03 03) 03 03) 03 03)
02 02) 02 02) 02 02)
03 02) 03 02) 03 02)
03 02) 03 02) 03 02)
04 03) 04 03) 04 03)
-1.4370048 (6.401786E -1.436426E (6.4019188 -1.4364028 (6401924E
-
-
-
+1.984058E (-2.249696E +1.983929E (-2.249736~ +1.983915E (-2.2497458 +4.1213278 (3.8891038 +4.1216318 (3.8886048 +4.1216288 (3.888565E
-
-
-2.273306E (-2.8273568 -2.272394E (-2.826639E -2.2723928 (-2.826613~
+6.327702E (6.8015558 +6*327264E (6.803787E +6*327223E (6.803575E
02 02) 02 02) 02 02)
03 03) 03 03) 03 03)
03 03) 03 03) 03 03)
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
used. Only thefirst 5 x 5 submatrix
w = 15.0, 0, = a/2, b = 1.0, k = 2, R = 3.0. The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
04 05) 04 05) 04 05)
02 02) 02 02) 02 02)
-
-2.021378E (-1.957353E -2.021449E (-1.957246E -2.021457E (-1.957234E
02 02) 02 02) 02 02)
-
+1*303077E (-1.009861E +1.303071E (-lGO9866E +1.303070E (-1.009872E
is shown below)
TABLE 1 Transmi$sion matrix T (There are$ve propagating modesfor the parameters
04 04) 04 04) 04 04)
05 04) 05 04) 05 04)
+4.88298OE (2.8368378 +4,883418E (2.835938E +4.883455E (2~835895E -
-
-
-
-8.69345OE (1.881421E -8.738582E (1.881236E -8.741247E (1.88123lE
+6.565637E (-lW7883E +6.561418E (-1 G47079E +6.561144E (--1.047054E
+6.786596E (-6.105665E +6.781811E (-6.107183E +6*781503E (-6.107251E
-
+l.O15965E (1.7999OlE +l.O15623E (1.802748E +l.O156OlE (1.80291SE
-
-
+6444452E (-2.244243E +6.438010E (-2.2507078 +6.4375948 (-2.25114OE
+8.15968OE (-1.247656E +8.1585108 (-1.250408E +8-l 58445E (-1.2505728
-5WO768E (-1.508793E -5WO774E (-1.512977E -5+300774E (-1.513187E
+4.769656E (2.770999E i4770084E (2,770122E +4,770120E (2.77008OE -
03 04) 03 04) 03 04)
04 05) 04 05) 04 05)
04 04) 04 04) 04 04)
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
-
-5443959E (-9822010E -5449999E (-9.821526E -5.450338E (-9.821465E
-
+6443038E (-9.300689E +6441785E (-9.26903 1E +6441669E (-9.267227E -
-4996537E (2.235231E -4.9965398 (2.235711E -4.996539E (2.235754E
+7.5482578 (-1.1541668 +7.547 174E (-1.156712E +7,547114E (-1.156864E -
-7.855391E (1.70005OE -7.896172E (1.699883E -7.898581E (1.699878E -
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
05 04) 05 04) 05 04)
-
-
-5.004256E (-7.1759OOE -5.004258E (-7.172510E -5.004258E (-7.172224E +8.900146E (-1.259086E +8.900764E (-1.258989E +8.90077OE (-1.258977E
-
-
+5.055236E (-1.760457E +5.050183E (-1.765527E +5+M9857E (-1.765867E +5.463521E (-7.8867328 +5*462459E (-7.859887E +5.462361E (-7.858357E
-
04 03) 04 03) 04 03)
01 04) 01 04) 01 04)
04 05) 04 05) 04 05)
04 05) 04 05) 04 05)
04 05) 04 05) 04 05)
used)
+5.03077OE (-8.029163E +5.027537E (--8.023005E +5.0273278 (-8.022809E
-5.001183E (-3,088994E -5.001183E (-3.089944E -5+M1183E (-3.08995OE
+5.990322E (-8.47439OE +5.990737E (-8.473738E +5.990742E (-8.473659E
-3.107061E (-5605771E -3.110508E (-5605495E -3.110702E (-5.60546OE
$5.363981 E (9.5029228 +5.362178E (9.517953E +5.362061E (9.518833E
+3,499956E (-3.148789E +3.497488E (-3.149572E f3.497330E (-3.1496078
-
-
-
-
01 04) 01 04) 01 04)
04 04) 04 04) 04 04)
04 04) 04 04) 04 04)
04 05) 04 05) 04 05)
___ - 04 - 04) - 04 - 04) - 04 - 04)
The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
04 04) 04 04) 04 04)
04 04 04 04) 04 04)
01 05) 01 05) 01 05)
-
-5~000759E (-8.849512E -5@00763E (-8.874411E -5.000763E (-8.875728E
TABLE 2 Matrix RI (see Table I .for values qf various parameters
06 05) 06 05) 06 05) 05 05) 05 05) 05 05)
-7.1443358 (-4.476370E -7.1388558 (-4.473493E -7.13868lE (-4.473307E -4.162767E (-4.627017E -4.163802E (-4.6237548 -4.163848E (46235458
05 02) 05 02) 05 02) 06 05) 06 05) 06 05)
-
-1.053117E (-3.489396E -1.0560378 (-3.489392E -1.056184E (-3.4893928 -8.7084678 (-5.695346E -8.727677E (-5.694529E -8.728827E (-5.694484E -1.566451E (-4.49814OE -1.570963E (-4.493644E -1.571265E (-4.493354E +1.2563078 (-7.0912978 +1.258294E (-7.088915E +1.258410E (-7.088759E
05 05) 05 05) 05 05)
06 05) 06 05) 06 05)
05 05) 05 05) 05 05)
+1*93412OE (-3.329156E +1.9335088 (-3.329454E +1.9334788 (-3.329479E -
06 05) 06 05) 06 05)
-
-7.0176OOE (-4.861431E -6.993713E (4.860486E -6.992352E (-4.860399E -7.4109498 (4.107602E -7.410584E (4.1121598 -7.4105388 (4.112415E
05 05) 05 05) 05 05)
05 06) 05 06) 05 06) 06 05) 06 05) 06 05) 05 02) 05 02) 05 02)
+1.282730E (5.927087E +1.282603E (5.957857E +1.2826OOE (5.9596748 -8.708468E (-5.6953468 -8.7276778 (-5.694529E -8.7288278 (-5.694484E +1.686537E (-3.7752378 +1.6868998 (-3.775235E +1.686932E (-3.775235E 05 02) 05 02) 05 02) 04 05) 04 05) 04 05)
-6.385103E (4445204E -6.382086E (4445203E -6.381832E (-4445203E -1.120332E (-7.919334E -1.120246E (-7.919883E -1.120236E (-7.919889E
-
06 05) 06 05) 06 05)
-
-7.0176OOE (4861431E -6.993713E (-4860486E -6.992352E (4.860399E
06 05) 06 05) 06 05) 06 05) 06 05) 06 05)
-
-7.1443358 (447637OE -7.138856E (4.4734938 -7.138681E (4473307E -1.566451E (449814oE -1.570963E (-4.493644E -1.571265E (-4.4933548
-
04 05) 04 05) 04 05) 4.083716E - 05 (-6.608544E - 02) -4.084972E - 05 (-6608543E - 02) 4084980E - 05 (-6608543E - 02)
-1.120332E (-7.919334E -1.1202468 (-7.919883E -1.120236E (-7.919889E
-4.162767E - 05 (4.627017E - 05) -4.163802E - 05 (4.6237558 - 05) -4.163848E - 05 (4.6235458 - 05) +1*256307E - 05 (-7.0912978 - 05) +1.258294E - 05 (-7.088915E - 05) +1.258410E - 05 (-7.0887598 - 05) -7*41095OE - 05 (4.107602E - 05) -7.410584E - 05 (4.112159E - 05) -7.410538E - 05 (4.112415E - 05)
The imaginary parts are enclosed in brackets. The first two numbers of each of the above matrix elements are for N = 10. The second and third pairs of numbers are for N = 20 and N = 25, respectively.
-
06 02) 06 02) 06 02) 05 05) 05 05) 05 05) 05 06) 05 06) 05 06)
-6.0334888 (-3408421E -6.050463E (-3408418E -6.0513618 (-3408418E +1.934120E (-3.329156E +1.933508E (-3.329454E +I .933478E (-3.329479E +I .28273OE (5.927087E +1.2826048 (5*957857E +1+2826OOE(5.959675E -
TABLE 3
Matrix R, (see Table 1 for values of various parameters used)
102
C. K. W. TAM TABLE 4
Eigenvalues (A:) w = 40x, /I = 1.O, k = 11, R = 1 (only thefirst 40 eigenvalues are tabulated below) N= 100
N=50
__1
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
-6.003698E -5.280253E -4.8755768 -4.572737E -4.3 10502E -4.079362E -3.869206E -3.676555E -3.497749E -3.330858E -3.174172E -3.026407E -2.886661E -2.753989E -2.62788OE -2.507594E -2.392832E -2.282982E -2.177865E -2.0769538 -3.980124E -1.886885E -1.797106E -1.710183E -1.625543E -1.541626E -1.459617E -1.406372E -1.34694OE -1.25248OE -1.148122E -1.0371OOE -9.202061E -7,978 184E -6.701889E -5.374953E -3.99876OE -2.574394E -1.102737E +4 154960E
+ 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 -t 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 03 + 03 + 03 + 03 -t- 03 -t 03 + 03 + 02
-6.003699E -5.280254E -4.875576E -4.572737E -4.3 10502E -4.079362E -3.869206E -3.676555E -3.497749E -3.330858E -3.174172E -3.026407E -2.886661E -2.753989E -2.627880E -2.507594E -2.39283 1E -2.282982E -2.177865E -2.076953E -1.980124E -1.886885E -1.797106E -1.710183E -1.625543E -1.541626E -1.459616E -14063728 -1.346941E -1.25248OE -1.148122E -1.037lOOE -9.20206OE -7.978 I80E -6.791885E -5.374946E -3.998753E -2.574386E -1.102736E +4.154727E
+ 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 + 04 -I- 04 -I- 04 + 04 f 04 + 04 + 04 + 04 + 04 + 04 + 04 + 03 + 03 + 03 + 03 + 03 + 03 + 03 + 02
of practical interest. With the above ordering of eigenvalues the relevant information about these modes are contained in the upper left-hand submatrices of R and T. Thus only these submatrices need be considered in actual computation. Tables 1, 2 and 3 represent some typical results of a routine computation. The numerical values of the parameters used are w = 15.0, p = 1.0, k = 2, f3, = 71/2 and R = 3.0. In this case there are five negative eigenvalues and hence five propagating modes. The first 5 x 5 submatrices of the matrices T, RI and Rz are given in the tables. Three sets of values with sizes of the matrices truncated to N = 10, 20 and 25 are exhibited. It can be seen that there is a rapid numerical convergence using the present computation procedure. Table 4 provides an example of numerical results for computations involving high frequency and high mode number incident wave. The parameters are w = 4071, ,9 = 1.O, k = 11 and R = 1
SOUND TRANSMISSION IN DUCT BENDS
103
25 -
w
Figure 3. Eigenvalues of propagating
modes. k = 2; fi= 1.0; R = 3.0.
(very sharp bend). There are thirty-nine negative eigenvalues and hence thirty-nine propagating modes (with the same k value). Given in the first column of this table are the first forty eigenvalues computed by using truncated matrices of 50 x 50. They are to be compared with the first forty eigenvalues obtained by using matrices of 100 x 100 shown in the second column. It is clear from such a comparison that the negative eigenvalues are accurate to at least six significant figures. This accuracy is sufficient for most applications. Figure 3 provides an example illustrating the dependence of Im(&), or the angular wave number, on the frequency, o. Only the eigenvalues of the propagating modes are shown. As frequency increases there are more and more propagating modes just as in the case of a straight duct. It is to be noted that even at moderately large values of o the angular wave numbers of the propagating modes are quite large. This implies (verified by actual computation) that the phases ofthe transmission and reflection coefficients will oscillate very rapidly between positive and negative values even with small changes in 19~or CD.Thus in considering the transmission of sound through a duct system good accuracy must be retained in joining the solutions together at the two ends of a curved bend so as to ensure correct matching of phases. 5. DISCUSSION
To make the above computation procedure practical it is necessary to truncate the matrices to a finite size of N x N. Naturally, the larger the N is the more accurate the computed solution will be. For some problems, however, only a reasonably accurate solution is all that is required. In these cases it is meaningful to ask what is the smallest N for which such a solution could be expected. To estimate the value of Nminit is worthwhile to note that the eigenvalues of nt are real. Further the negative eigenvalues give rise to propagating modes while the positive eigenvalues provide the cutoff modes. For a bend with a not too small turning angle dE, the cut-off modes play only a very insignificant role in the transmission of sound. They are rapidly attenuated spatially. This observation indicates, therefore, that Nmin is an integer somewhat larger than the total number of negative eigenvalues. This is confirmed computationally. An illustration of this is given by the very accurate results of the first column of Table 4 where N = 50 (39 negative eigenvalues) is used.
104
C. K. W.
TAM
One interesting
phenomenon observed in our numerical experiments testing the present procedure is that the transmission matrix T is symmetric (see Table 1). From the second of equations (17) it follows that T is also symmetric and that the index j of the incident wave appears only in the diagonal elements of T. The physical implication is that a reciprocal relationship between incident and transmitted wave modes exists. That is to say that the transmission coefficient for the mth modes due to an incident wave of the nth mode is the same as the transmission coefficient of the nth mode due to an incident wave of the mth mode. This
interesting property does not seem to have been noticed before. In the work reported here the computation procedure has been developed primarily for rectangular ducts with rigid walls. However, it is straightforward to see that with only simple modifications the same procedure could be used for ducts with soft walls. Moreover, by using a toroidal co-ordinate system instead of a cylindrical co-ordinate system the present method is applicable to the study of sound transmission in circular bends of circular ducts. In this case the wave equation is non-separable so that the classical method of separation of variables cannot be used. It is envisaged that computation results involving soft walled rectangular ducts and circular ducts will be reported in a subsequent paper. ACKNOWLEDGMENT This work was supported in part by the National Science Foundation under Grant GK35790 and the National Aeronautical and Space Administration Langley Research Center Grant NSG 1021. The author wishes to thank Barbara Schreur for her assistance in computer programming. REFERENCES 1. P. E. KRASNUSHKIN 1945 Uch. Zap. Mosk. Gos. Univ. No. 75, Bk. 2, Pt. 2, 9-27.
On waves in curved tubes. 2. A. G. SEVESHNIKOV1958 Radiotekhnika i Eiektrouika 3, 641-648. Waves in bent tubes. 3. F. E. GRIGOR’YAN 1969 Akustische Zeitschrift 14, 376-384. (English translation: Soviet Physics Acoustics 14, 315-321, 1969). Theory of sound wave propagation in curvilinear wave guides. 4. W. C. OSBORNE1968 National College for Heating, Ventilating, Refrigeration and Fan Engineering Technical Memo. No. 5. The propagation of sound waves round a bend of rectangular section. 5. W. ROSTAFINSKI1972 Journalof the AcousticalSociety ofAmerica52,141 I-1420. On propagation of long waves in curved ducts. 6. W. ROSTAFINSKI1974 Journal of the Acoustical Society of America 56,l l-15. Analysis of propagation waves of acoustic frequencies in curved ducts. 7. A. CUMMINGS1974 Journal of Sound and Vibration 35, 451-477. Sound transmission in curved duct bends. 8. M. K. MYERS and P. MUNGUR 1975 American Institute of Aeronautics and Astronautics Paper 75-497. Sound propagation in curved ducts. 9. B. A. FINLAY~~N1972 The Method of Weighted Residuals and Variational Principles. New York: Academic Press. 10. W. EVERSMAN,E. L. COOK and R. J. BECKEMEYER1975 Journal of Sound and Vibration 38, 105-123. A method of weighted residuals for the investigation of sound transmission in nonuniform ducts without how. 11. J. H. WILKINSON 1965 The Algebraic Eigenvalue Problem. The Oxford University Press.