Transfer matrix of a truncated cone with viscothermal losses: application of the WKB method

Summary The propagation in tubes with varying cross section and wall visco-thermal eﬀects is a classical problem in musical acoustics. To treat this aspect, the ﬁrst method was the division in a large number of short cylinders. The division in short conical frustums with wall eﬀects independent of the radius is better, but remains time consuming for narrow tubes and low frequencies. The use of the WKB method for the transfer matrix of a truncated cone without any division is investigated. In the frequency domain, the equations due to Zwikker and Kosten are used to deﬁne a reference result for a simpliﬁed bassoon by considering a division in small conical frustums. Then expressions of the transfer matrix at the WKB zeroth and the second orders are derived. The WKB second order is good at higher frequencies. At low frequencies, the errors are not negligible, and the WKB zeroth order seems to be better. This is due to a slow convergence of the WKB expansion for the particular case: the zeroth order can be kept if the length of the missing cone is large compared to the wavelength. Finally, a simpliﬁed version seems to be a satisfactory compromise.


Introduction
The calculation of the transfer matrix of wind instrument resonators is an old problem. The most general model is one-dimensional, based on the horn equation, written for either plane or spherical waves. A classical difficulty is the effect of boundary layers in tubes with variable cross section area, which depends on the radius. For this reason, no analytic expression for the transfer matrix of a truncated cone with visco-thermal losses was established yet. The general idea is to divide the resonator into frustums of cylinders (see [1]). In order to limit the time consumption, another method uses conical frustums with boundary layer effects independent of the radius [2,3,4]. The aim of the present paper is to use the well known WKB (Wentzel-Kramers-Brillouin) method, in order to reduce the computing time by limiting the division of a truncated cone to one segment. For this purpose, it is necessary to state the wave equation with visco-thermal effects without term of first-order (space) derivative. Reducing the computing time is not useful for a single input impedance computation, but this becomes useful for applications such as optimization processes, when numerous input impedance computations are requested.
All calculations are done in the frequency domain. In Section 2, the classical derivation of the transfer matrix of a lossless, truncated cone is briefly reminded, for spherical and plane waves. In Section 3 the result of the Zwikker and Kosten (ZK) theory [5] is used. The second order with respect to the Stokes number is considered. The Stokes number is the ratio of the radius to the boundary layer thickness, which is inversely proportional to the square root of the frequency.
In Section 4 reference numerical results are sought by studying their convergence when dividing the cone with an increasing number of conical frustums. Because losses (and dispersion) are particularly strong in instruments with low frequencies and narrow radii, the input impedance of a simplified bassoon is the main example considered.
Then, in Section 5, the WKB method is used in order to find approximate solutions of the Helmholtz equation with losses. In Section 6, these solutions lead to an original expression of the transfer matrix, which allows the number of segments to be reduced to one. In Section 7 the different formulas are compared and discussed. Section 8 proposes simplified, approximate formulas for both the spherical and the plane wave approximation, and Section 9 presents a conclusion.
2 Helmholtz equation and transfer matrix for truncated cones 2

.1 Geometry
Consider a truncated cone of length L, input radius R 1 , output radius R 2 . If a spherical surface (wavefront) is assumed, the curvilinear abscissa r is preferred. The radius is R = r sin ϑ, where ϑ is the half-angle at the apex. The length L = r 2 − r 1 is equal to (R 2 − R 1 )/ sin ϑ. The area of the spherical wavefront is Σ = 2πr If a planar wavefront is assumed, the longitudinal coordinate denoted x is used, with R = x tan ϑ. The length = x 2 − x 1 is equal to L/ cos ϑ. The area of the planar wavefront is simply S = πR 2 . For small angles, the two expressions of the area are equivalent, at the second order of ϑ.

Spherical waves
In the frequency domain, k = ω/c is the wavenumber (c is the speed of sound, ω the angular frequency). For spherical wavefront, the Helmholtz equation for the acoustic pressure P (r) is exactly derived from the 3D Helmholtz (lossless) equation. It is written as: and the general solution is: Using the Euler equation, ∂ r P = −jkρcV , where ρ is the air density, the particle velocity V is given by The flow rate U = ΣV is deduced. Two intermediate variables are defined: with the following relationships: The product of the 3 matrices lead to the following expression of the transfer matrix: The value of C is obtained by using reciprocity. It can be rewritten in order to use directly for the two kinds of cones (either diverging or converging): The second factor tends to unity for cylinders (R/r tends to zero). In what follows, it is replaced by this limits in the expression of the transfer matrices.

Plane waves
If plane wavefronts are assumed, the previous derivation remains valid if r is replaced by x, L is replaced by , and Σ is replaced by S. This yields: This expression was given in [8], and in another form, in [9].

ZK theory and Helmholtz equation with losses
3.1 Equations of the ZK theory for cylindrical tubes; extension to plane waves in conical tubes In this paper, we use the Zwikker and Kosten theory [5], written in the form of a telegraphist equation, such as: For a cylinder, the parameters by unit length are Z v and Y t the series impedance and shunt admittance, respectively. They are constant. The well known expression of the transfer matrix is the following: The wavenumber K and the characteristic impedance Z c are: The theory is based upon the (excellent) approximation that the pressure is planar, independent of the abscissa x . The parameters are: with: (22) γ h is the ratio of the specific heats. k v and k t are the viscous and thermal wavenumbers, respectively: µ is the viscosity, κ the thermal conductivity, and C p is the specific heat at constant pressure. The Stokes number is denoted S t = |k v | R . At the second order of S t , the asymptotic expression of Eq. (22) to the viscous parameter is: We need also the value of y t , or that of the wavenumber, which is denoted K = kq. The asymptotic expansion is well known (see e.g. [10,8]). At the same order of S t where For a cone, we assume that at each abscissa of a cone, the pressure is plane and the ZK theory is valid. The parameters Z v and Y t depend of the abscissa x, and the transmission line is non-uniform. The theory is assumed to be valid for tubes of variable cross section S(x), with transmission lines parameters Z v and Y t . It is possible to write: where, because R = x tan ϑ, In the literature on musical acoustics, the Euler equation is often kept unchanged (i.e., z v = 1, see e.g. [6,7]). Another approximation if often proposed: the characteristic impedance is assumed to be constant, i.e. z v = y t = q [8,9]. Moreover the expressions are in general limited to the first order of the Stokes number.

Differential equation for the pressure in a tube of variable cross section
We are searching for a generalization of Eqs.(14 to 16). For a tube with variable cross section, if planar wavefronts are assumed, the pressure equation can be written, from Eqs.(17, 18), as: where the prime symbol indicates the derivative with respect to the abscissa x. By changing the variable, this Sturm-Liouville equation can be rewritten into a canonical form (without the first derivative of P ), as follows: Using the second derivative of Eq. (31), it is obtained: Together with Eq. (32), this leads to : the coefficient η is expressed as follows: For a cone, R = 0, R = x tan ϑ, and R /R = 1/x . Thus: As a result, for a cone, This simple result is a particularity of the axisymmetrical geometry. Furthermore, at the first order of α/x, the vanishing of η is valid for any shape of circular tube.

Differential equation for plane waves in a cone
Finally, for a cone, the differential equation to be solved is, at the second order of α/x: where and q is given by Eq. (26). This original equation is extremely simple. In the literature concerning the WKB solution of the Schrödinger equation the term k −2 is the analogous to the square of the Planck constant .
4 Computation of a reference solution.

Numerical, "exact" results with division into conical frustums
In order to evaluate the errors made using the WKB solution, a numerical, "exact" solution is sought. The only geometry for which there is an analytical solution for the Helmholtz equation with losses is the cylinder. It is possible to approach the cone by cylindrical frustums of very short length and constant radius R n . Eqs. (19) give the transfer matrix of a cylindrical tube of length . However the convergence with the number of frustums is more rapid if conical frustums with constant losses are used. For that purpose the cone is divided in several short conical frustums with a small radius variation. The losses are approximated by the losses of a cylinder with an equivalent radius R eq . In order to correctly estimate the magnitude of the impedance peaks, Ref. [8] suggests, for the frustum between x 1 and x 2 , to choose this radius as follows: where = x 2 − x 1 is the length of the considered frustum. Losses are computed by using Eqs. (25 and 26) with the equivalent radius R eq . The transfer matrix of a conical frustums is given by Eqs.(14 to 16), K and Z c being calculated for the the equivalent radius. The example of a simplified bassoon (R 1 = 0.002 m; R 2 = 0.02 m; L = 2.43 m) is chosen for the low values of the Stokes number at low frequencies: for the first resonance frequency, it is close to 10, and its inverse is 0.1. The input impedance is computed by multiplying the transfer matrices of each slice and by assuming a zero impedance at the wide end (non-radiating open pipe) (Fig.2). The computation is carried out for the frequency range [20, 10 4 ]Hz, with a logarithmic step of 1 cent. The cent is a musical logarithmic scale (100 cents = 1 semitone) defined as dif (cents) = 1200 * log 2 (f 2 /f 1 ). In order to determine the frustum length necessary to converge to the "exact" solution, the evolution of the norm of the difference to the finest slicing is observed (Fig.3) The norm converges for slices of 0.1mm corresponding to about 2 · 10 4 frustums for the considered tube. With this slicing the variation of the radius for one frustum is about 7 · 10 −7 m. Throughout this paper, the reference impedance Z ref corresponds to that computed with the finest slicing (10 5 conical frustums of 2.43 · 10 −5 m).
Conversely, Fig.3 shows that the model with cylindrical frustums seems to converge also toward the reference impedance (Fig.3), but the convergence is slow and a too large number of cylinder would be necessary to converge (at least 10 10 ). Even if the convergence is not reached, the norm of the difference between this theoretical results and the reference impedance for the finest slicing is very small (10 −4 ). The norm of the impedance difference is not the more convenient way to observe the difference between impedance curves. The difference between the resonance peak characteristics is a more significant observation, especially for musical instruments. The resonance frequencies f (i) and magnitudes a (i) of the impedance peaks are computed with the finest cylindrical slicing and compared to those of the reference impedance in Figure 4. They are estimated by applying a second order polynomial fit on the modulus in decibel over the three samples around the peaks to be more precise than the frequency step. The resonance frequency deviations, given in cents, are under 0.01 cents for all peaks which is negligible for musical application (it is generally assumed that the human ears can not detect frequency difference smaller than 3 cents) (Fig.4, a). The peaks magnitude is very well estimated with a difference within 1 · 10 −3 dB (Fig.4, b).

Computation using asymptotic expressions with respect to the Stokes number
The impedance is also computed for conical frustums with the asymptotic expression of z v and y t at the first and second order of the Stokes number (sec. 3.3). The norm of the difference with the reference impedance decreases with the slice length to reach 10 −2 for the first order and 3 · 10 −4 for the second order (Fig.3). These values seem high but it is important to notice that the norm of the error is very sensitive and particularly to a small frequency decay. Figs. 3 and 4 show that the resonance frequency deviations are within 0.1 cents for all approximations (Fig.4, a). This is acceptable in a musical context. The magnitudes also are well estimated by the different approximations, even if the first order approximation has a significantly lower accuracy (the deviation is 30 times higher in dB) than the second order approximation (Fig.4, b).
As a conclusion, the asymptotic expression of the visco-thermal losses at the second order of the Stokes number is sufficient to achieve a satisfactory accuracy.

WKB solution of the second order equation
We present the classical derivation of the WKB method (see Ref.citevoros ). The solution ψ of Eq. ( 41) is sought in the following form: ψ = g e jk udr (49) and its derivative is equal to: The unknowns are the functions g and u. u is dimensionless. Calculating the second derivative ψ", the following equation is obtained: Because they are two unknowns functions, the vanishing of the bracket can be chosen, This leads to the following relationship: where Λ is a constant. The first equation yields: If u is a solution, −u is also a solution. Therefore ψ is the superposition of two solutions, as follows: where ψ ± are two constants. u is the unique unknown function. It is solution of the following equation, which is derived from Eqs. (52) and (53): For infinite k, a trivial solution is u = q. For large k, the solution is sought in the form of an asymptotic expansion: Simplifying by the factor 1/k 2 , the term v can be obtained: therefore, at the second order of 1/k : Another useful expression for the function v is: At the fourth order of 1/k, w is obtained: The expression of the fourth order coefficient leads to complicated results and it is not studied in the present paper.
6 Transfer matrix for the WKB solution 6

.1 Expression of the transfer matrix
The derivation of the transfer matrix is similar to that for the non dissipative case (Section 2) Using Eqs. (31) and (54), the relationship between the pressure P andP is found to be: The variableÛ is defined as:Û For the variablesP andÛ , the transfer matrix is written as : This integration can be interpreted as an average of the visco-thermal effects over the conical frustum. The development of the expression of u makes appearing the integration of 1/R (Eqs. (28, 56)), which is coherent with the equivalent radius R eq chosen in Eq. 43. The relationship between (P, U ) and (P ,Û ) is given by Eq. (17): Thus, with Eq. (63): The transition matrix at each extremity of the truncated cone is derived as follows: with The determinant of the matrix is a i d i = jk and is independent of the extremity (1 or 2). The inverse matrix is: The product of the 3 matrices between abscissae 1 and 2 gives the final transfer matrix: The determinants of the transition matrix and its inverse matrix are inverse. Therefore the determinant of the transfer matrix is unity, and the coefficient C can be derived from the other coefficients. The final result is: As a summary, the calculation of the transfer matrix requires the values of the definite integral of the function u, of the function u itself and of its derivative. Consequently, according to Eq. (56), we need the equivalent expressions for the functions q and v, their integral and their derivative.

Expansions of q and v with respect to the Stokes number
The starting point is the expansion of the quantity q, as written in Eq. (26), at the second order of the Stokes number. The indefinite integral I q = k x2 x1 qdx and the derivative q are given by: Notice that x is a variable with dimension, thus the neperian logarithm ln(x) has no sense, but, in the transfer matrix, the quantity ln(x 2 /x 1 ) intervenes in X. The function v could be also expanded with respect to α/x, but for the integral I v it is simpler to use directly Eq. (59). The integral of the second term of the expression (59) is complicated, but a numerical evaluation is found to be very small, therefore In order to calculate v , we need the calculation of the second and third derivatives of q (Eq. (58): Finally, we need to successively compute the following expressions: 7 Comparison of the zeroth and second orders of the WKB solution

Discussion
We wish to compare the "exact" results with the WKB solutions at zeroth and second orders in k −1 . The zeroth order is obtained by writting u = q. The comparison is done with the reference result obtained in Sect.

4.1.
We first notice that the obtained transfer matrix given by Eqs. (74 to 76) is invariant by slicing: the product of the matrix of two successive conical frustums with the same angle equals the transfer matrix of the total cone. This invariance can be shown by a rather long calculation for both the orders of approximation. It is not difficult, and is not developed here.
The analytical limit of validity of the WKB solution can be given by the comparison between the terms due to the visco-thermal effects and their variations with the radius. At the first order of the Stokes number and at the second order of the WKB solution, the function u is written as: The second term in the bracket can be very high at low frequencies for the smaller radius of the cone. For the first resonance of the bassoon, it is close to 5, and it appears that the convergence to a correct limit value at higher WKB order is problematic. Similar behaviour can be found for the integral and the derivative: and Using Eq. (82) the quantity X (Eq. (67) can be calculated. The result is: Therefore, an approximate condition of validity can be estimated: kx 1 > 1 . For the first peak of a divergent cone, k ≈ π/x 2 , thus the condition can be written as: This is not satisfied for the first peak of a bassoon. However, above the third or fourth peak, this becomes acceptable.

Numerical comparison
The impedances are computed with the zeroth order (u = q) and 2nd order (u = q + v/k 2 )of the WKB transfer matrix for a truncated cone. They are compared to the reference impedance (average ZK losses for 10 5 conical frustums) in Figure 5,a. As expected, the WKB approximations are better at higher frequencies. It is particularly true for the second order approximation, for which the first peak is underestimated and the second peak largely overestimated. Both correspond to frequencies for which kx 1 < 1. However, the zeroth order approximation seems better at this peaks. This observation can be quantified by applying the error norm of Eq.(48) of different frequency ranges. If the norm is computed on the entire range, the second order result has poorer accuracy than the zeroth order one (the error norm are respectively 0.3 and 0.1). Conversely, if the range is limited to higher frequencies (kx 1 > 2), the second order is better than the zeroth order (respectively 0.004 and 0.01).
These observations appear more clearly on the modal parameters of the impedance peaks, which are represented versus both the frequency and the normalized wavenumber kx 1 in Figure 6. At low frequencies (kx 1 < 1.5, Fig.6,a) the zeroth order is better. It induces a deviation about −30cents and −4dB on the first peak, while the second order induces a deviation about 140 cents and −14 dB. As expected, at higher frequencies (kx 1 > 2, Fig.6, b), the second order is slightly better, even if both orders show acceptable deviations to the reference impedance (< 1 cents and < 0.5dB).
Consequently the gain of the second order of WKB on the zeroth order is questionable for this kind of resonators. The zeroth order approximation seems much more simpler and more appropriate for musical applications, when the first peak is predominant in the sound production. It can be useful to search for a further simplification of the zeroth order solution. We start from Eqs. (74 to 76) with q 1 = u 1 and q 2 = u 2 . First of all, we notice the following relationship: Z v /q = jkZ c . The transfer matrix can be written as follows: At the second order of the Stokes number, the characteristic impedance is written as: Therefore, by comparing with the expression of the wavenumber K = kq (see Eq. (28)), the visco-thermal effects intervene significantly more in K than in √ Z c (1.044 and 0.18, respectively). For this reason, we can ignore the influence of these effects in Z c (This is due to the difference in sign of viscous and thermal effects. Moreover, in the transfer matrix expression, the characteristic impedance intervenes through its variation in the cone, between x 1 and x 2 ). As a consequence, the following simplified formula can be useful: (91) Figure 7 shows that the simplification leads to satisfactory results, when compared to the complete zeroth order formula. Figure 7: Difference between the resonance parameters of the WKB simulations at the zeroth order (complete and simplified) to those of the reference impedance: top: frequencies and bottom: magnitudes. In both simulations, the second order approximation is used for the visco-thermal effects.
This approximation differs from other approximations published by some authors. The result of Ref. [11] adds the approximation that q 1 = q 2 = 1, and does not justify the starting point. Other approaches were given for the calculation of the impedance peaks [8,13] or of the resonance frequencies [14], the results being also consistent with the equation (67), which is based upon the calculation of a mean value of the inverse of the radius.
We emphasize that this paper considers an extreme case, the bassoon. Fig. 8 shows that for a wider instrument, such as a simplified baritone saxophone (R 1 = 6.75mm, R 2 = 60mm, = 2.38m), the final results are much better for the first peaks. Furthermore the convergence study to the reference result justifies a classical approximation of visco-thermal effects in cones: the second order of the asymptotic expansion with respect to the Stokes number does not significantly improve the accuracy of the computation. This was not clear in [13].

Spherical waves
The previous analysis assumes plane waves. Nevertheless, for cones with wide apex angle, it is known that spherical waves are more accurate. This is especially important for the division of bells into truncated cones [7]. For such cones, visco-thermal effects are weak, and the WKB solution of zeroth order is satisfactory. It is s not necessary to repeat the complete analysis. We propose the following formula, which is a modification of Eqs. (8 to 10): with and I q = r + γαln(r) − δ α 2 r . (94)

Conclusion
To our knowledge the propagation equation in cones with visco-thermal effects was not derived in previous papers. In the frequency domain, at the second order of the asymptotic expansion with respect to the Stokes number, the canonical form (without term with the derivative of the first order) is extremely simple. With this starting point, the WKB method leads to the possibility to compute the transfer matrix of a truncated cone without division of its length. The result of the zeroth order is satisfactory under the condition that the length of the missing cone x 1 is large compared to the wavelength. Under this condition, the WKB second order is even better. However, the good accuracy of these approximations at higher frequencies is not actually necessary in practice. Conversely, for the first resonances, the accuracy is not excellent. Paradoxically the second order seems less accurate than the first one. This is due to the slow convergence of the series expansion of the WKB function denoted u in this paper. It should be necessary to extend the WKB method to further orders, but this should lead to complicated expressions. Therefore, we propose to limit the expansion to the zeroth order WKB formula, which is well known. Moreover, from an analysis of the dependence of the characteristic impedance with respect to the Stokes number, a significant simplification of the zeroth order is obtained. It slightly differs from formulas found in the literature, but a numerical analysis of the analytic formulas shows that the errors of the various formulas are of the same order of magnitude.
The values of the transfer matrix coefficients can be empirically improved at low frequencies. This is true in particular for the examined case of the input impedance of the simplified bassoon: ignoring the visco-thermal effects in the coefficients q 1 and q 2 (i.e., writing q 1 = q 2 = 1) diminishes the errors by a factor 2 for the first two peaks. However, this do not seem to be general for all coefficients of the matrix. We prefer to keep the formulas analytically derived with a clear expression of their origin. Formulas (92) are a satisfactory compromise for conical instruments.
Finally, these unified formulas for cones with either narrow and wide apex angles can be appropriate for impedance computation software. This is particularly useful for the application to instruments with bells, when the plane wave approximation is not applicable. An extension of the present work could be taken for other shapes of horns, as it was done in [7], but the corresponding work should be heavy.