Issue 
Acta Acust.
Volume 4, Number 2, 2020



Article Number  7  
Number of page(s)  11  
Section  Musical Acoustics  
DOI  https://doi.org/10.1051/aacus/2020005  
Published online  12 May 2020 
Technical & Applied Article
Transfer matrix of a truncated cone with viscothermal losses: application of the WKB method
^{1}
Magique 3D Team, Inria Bordeaux Sud Ouest, 200 Avenue de la Vieille Tour, 33405 Talence Cedex, France
^{2}
Aix Marseille University, CNRS, Centrale Marseille, LMA UMR 7031, 13013 Marseille, France
^{*} Corresponding author: augustin.ernoult@gmail.com
Received:
20
January
2020
Accepted:
23
April
2020
The propagation in tubes with varying cross section and wall viscothermal effects is a classical problem in musical acoustics. To treat this aspect, the first method is the division in a large number of short cylinders. The division in short conical frustums with uniform averaged wall effects 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 define a reference result for a simplified 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 simplified version seems to be a satisfactory compromise.
Key words: Conical tubes / Viscothermal effects / WKB method
© A. Ernoult and J. Kergomard, Published by EDP Sciences, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The calculation of the transfer matrix of wind instrument resonators is an old problem. The most general model is onedimensional, 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 viscothermal 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 equal to those of a cylinder of similar length and equivalent radius [2–4]. The aim of the present paper is to use the well known Wentzel–Kramers–Brillouin (WKB) 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 viscothermal effects without term of firstorder (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 inverse of 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} (Fig. 1). If a spherical surface (wavefront) is assumed, the curvilinear abscissa r is preferred. The radius is R = r sin ϑ, where ϑ is the halfangle at the apex. The length L = r_{2} − r_{1} is equal to (R_{2} − R_{1})/sin $\vartheta $. The area of the spherical wavefront is $\mathrm{\Sigma}=2\pi {r}^{2}\left[1\mathrm{cos}\vartheta \right]=2\pi {R}^{2}/\left[1+\mathrm{cos}\vartheta \right]$.
Figure 1 Sketch of a conical tube. $L={r}_{2}{r}_{1}=\mathcal{l}/\mathrm{cos}\vartheta .$ 
If a planar wavefront is assumed, the longitudinal coordinate denoted x is used, with R = x tan $\vartheta $. The length $\mathcal{l}$ = x_{2} − x_{1} is equal to L cos $\vartheta $. 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 $\vartheta $ (Fig. 1).
2.2 Spherical waves
In the frequency domain, k = ω/c is the wavenumber (c is the speed of sound, ω is 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:
$$\frac{{\mathrm{d}}^{2}\left(\mathrm{rP}\right)}{\mathrm{d}{r}^{2}}+{k}^{2}\left(\mathrm{rP}\right)=0,$$(1)and the general solution is:
$$\mathrm{rP}={Q}^{+}{e}^{\mathrm{jkr}}+{Q}^{}{e}^{\mathrm{jkr}}.$$(2)
Using the Euler equation, ${\mathrm{\partial}}_{r}P=\mathrm{jk\rho}\mathrm{cV}$, where ρ is the air density, the particle velocity V is given by,
$$\mathrm{V\rho}\mathrm{cr}={Q}^{+}{e}^{\mathrm{jkr}}{Q}^{}{e}^{\mathrm{jkr}}+\frac{P}{\mathrm{jk}}.$$(3)
The flow rate U = ΣV is deduced. Two intermediate variables are defined:
$$\begin{array}{l}\widehat{P}={Q}^{+}{e}^{\mathrm{jkr}}+{Q}^{}{e}^{\mathrm{jkr}}\\ \widehat{U}={Q}^{+}{e}^{\mathrm{jkr}}{Q}^{}{e}^{\mathrm{jkr}},\end{array}$$with the following relationships:
$${\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{1}=\left(\begin{array}{cc}\mathrm{cos}\mathrm{kL}& j\mathrm{sin}\mathrm{kL}\\ j\mathrm{sin}\mathrm{kL}& \mathrm{cos}\mathrm{kL}\end{array}\right){\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{2},$$(4)
$${\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{\mathrm{1,2}}=\left(\begin{array}{ll}r& 0\\ \frac{1}{\mathrm{jk}}& r\frac{\rho c}{\mathrm{\Sigma}}\end{array}\right){\left(\begin{array}{l}P\\ U\end{array}\right)}_{\mathrm{1,2}},$$(5)
$${\left(\begin{array}{l}P\\ U\end{array}\right)}_{\mathrm{1,2}}=\left(\begin{array}{ll}\frac{1}{r}& 0\\ \frac{\mathrm{\Sigma}}{\rho c}\frac{1}{\mathrm{jk}{r}^{2}}& \frac{\mathrm{\Sigma}}{\rho c}\frac{1}{r}\end{array}\right){\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{\mathrm{1,2}}.$$(6)
The product of the three matrices lead to the following expression of the transfer matrix:
$${\left(\begin{array}{l}P\\ U\end{array}\right)}_{1}=\left(\begin{array}{ll}A& B\\ C& D\end{array}\right){\left(\begin{array}{l}P\\ U\end{array}\right)}_{2};$$(7)
$$A=\frac{{r}_{2}}{{r}_{1}}\mathrm{cos}\left(\mathrm{kL}\right)\frac{\mathrm{sin}\left(\mathrm{kL}\right)}{k{r}_{1}};$$(8)
$$D=\frac{{r}_{1}}{{r}_{2}}\mathrm{cos}\left(\mathrm{kL}\right)+\frac{\mathrm{sin}\left(\mathrm{kL}\right)}{k{r}_{2}};$$(9)
$$B=\frac{\varrho c}{{\mathrm{\Sigma}}_{2}}\frac{{r}_{2}}{{r}_{1}}j\mathrm{sin}\left(\mathrm{kL}\right);$$(10)
The value of C is obtained by using reciprocity. It can be rewritten in order to be used directly for the two kinds of cones (either diverging or converging):
$$A=\frac{{R}_{2}}{{R}_{1}}\mathrm{cos}\left(\mathrm{kL}\right)\frac{\mathrm{sin}\left(\mathrm{kL}\right)}{\mathrm{kL}}\frac{{R}_{2}{R}_{1}}{{R}_{1}};$$(12)
$$D=\frac{{R}_{1}}{{R}_{2}}\mathrm{cos}\left(\mathrm{kL}\right)+\frac{\mathrm{sin}\left(\mathrm{kL}\right)}{\mathrm{kL}}\frac{{R}_{2}{R}_{1}}{{R}_{2}};$$(13)
$$B=\frac{\rho c}{\pi {R}_{1}{R}_{2}}\frac{1+\mathrm{cos}\vartheta}{2}j\mathrm{sin}\left(\mathrm{kL}\right).$$(14)
The second factor of B in equation (14) 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.
2.3 Plane waves
If plane wavefronts are assumed, the previous derivation remains valid if r is replaced by x, L is replaced by $\mathcal{l}$, and Σ is replaced by S. This yields:
$$A=\frac{{R}_{2}}{{R}_{1}}\mathrm{cos}\left(k\mathcal{l}\right)\frac{\mathrm{sin}\left(k\mathcal{l}\right)}{k\mathcal{l}}\frac{{R}_{2}{R}_{1}}{{R}_{1}};$$(15)
$$D=\frac{{R}_{1}}{{R}_{2}}\mathrm{cos}\left(k\mathcal{l}\right)+\frac{\mathrm{sin}\left(k\mathcal{l}\right)}{k\mathcal{l}}\frac{{R}_{2}{R}_{1}}{{R}_{2}};$$(16)
$$B=\frac{\rho c}{\pi {R}_{1}{R}_{2}}j\mathrm{sin}\left(k\mathcal{l}\right).$$(17)
This expression was given in [6], and in another form, in [7].
3 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:
$$\mathrm{d}P/\mathrm{d}x={Z}_{v}\left(x\right)U;$$(18)
$$\mathrm{d}U/\mathrm{d}x={Y}_{t}\left(x\right)P.$$(19)
The parameters by unit length are Z_{v} and Y_{t}, the series impedance and shunt admittance, respectively.
The cylinder case, for which they are constant, is first remembered. The well known expression of the transfer matrix is the following:
$$\begin{array}{l}A=D=\mathrm{cos}K\mathcal{l}\\ B=j{Z}_{c}\mathrm{sin}K\mathcal{l},\end{array}$$(20)where the wavenumber K and the characteristic impedance Z_{c} are given by:
$$K=j({Z}_{v}{Y}_{t}{)}^{1/2}\hspace{1em}\mathrm{and}\hspace{1em}{Z}_{c}=({Z}_{v}/{Y}_{t}{)}^{1/2}.$$(21)
The theory is based upon the approximation that the wavefront is planar, independent of the abscissa x. The parameters are:
$$\begin{array}{l}{Z}_{v}=\mathrm{j\omega}\rho {z}_{v}/S\\ {Y}_{t}=\mathrm{j\omega}{y}_{t}S/\left(\rho {c}^{2}\right),\end{array}$$(22)with:
$${z}_{v}={\left[1\frac{2}{{k}_{v}R}\frac{{J}_{1}\left({k}_{v}R\right)}{{J}_{0}\left({k}_{v}R\right)}\right]}^{1};$$(23)
$${y}_{t}=1+\left({\gamma}_{h}1\right)\frac{2}{{k}_{v}R}\frac{{J}_{1}\left({k}_{t}R\right)}{{J}_{0}\left({k}_{t}R\right)},$$(24)γ_{h} is the ratio of the specific heats. k_{v} and k_{t} are the viscous and thermal wavenumbers, respectively:
$${k}_{v}=\sqrt{\frac{\mathrm{j\omega}\rho}{\mu}};\hspace{1em}{k}_{t}=\sqrt{\frac{\mathrm{j\omega}\rho {C}_{p}}{\kappa}},$$(25)μ is the viscosity, κ the thermal conductivity, and C_{p} is the specific heat at constant pressure. This give $K=k\sqrt{{z}_{v}{y}_{t}}$ and ${Z}_{c}=\rho c/S\sqrt{{z}_{v}/{y}_{t}}$. The Stokes number is $\mathrm{St}=\left{k}_{v}\rightR$. For high Stokes number (wide duct and/or high frequency), at the second order of 1/St, the asymptotic expression of equation (23) to the viscous parameter is:
$${z}_{v}=1+2\sqrt{j}S{t}^{1}3\mathrm{jS}{t}^{2}.$$(26)
We need also the value of y_{t}, or that of q = K/k. The asymptotic expansion is well known (see e.g., [6, 8]). At the same order of St,
$$q=\sqrt{{z}_{v}{y}_{t}}=1+\gamma S{t}^{1}+\delta S{t}^{2},$$(27)where
$$\gamma =1.044\sqrt{2j};\hspace{1em}\delta =1.08j.$$(28)
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 nonuniform. 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, for high Stokes number:
$$\begin{array}{l}{z}_{v}=1+2\sqrt{j}\frac{\alpha}{x}3j\frac{{\alpha}^{2}}{{x}^{2}},\\ q=\sqrt{{z}_{v}{y}_{t}}=1+\gamma \frac{\alpha}{x}+\delta \frac{{\alpha}^{2}}{{x}^{2}},\end{array}$$(29)where, because R = x tan $\vartheta $, and so $\mathrm{St}=x\mathrm{tan}\vartheta \left{k}_{v}\right$,
$$\alpha =\frac{x}{\mathrm{St}}=\frac{1}{\mathrm{tan}\vartheta}\sqrt{\frac{\mu}{\omega \rho}}.$$(30)
In the literature on musical acoustics, the Euler equation is often kept unchanged (i.e., z_{v} = 1, see e.g., [9, 10]). Another approximation, the Heaviside condition, if often proposed: the characteristic impedance is assumed to be constant, i.e. z_{v} = y_{t} = q [6, 7]. Moreover the expressions are in general limited to the first order of the inverse of the Stokes number.
3.2 Differential equation for the pressure in a tube of variable cross section
We are searching for a generalization of equations (15)–(17). For a tube with variable cross section, if planar wavefronts are assumed, the pressure equation can be written, from equations (18) and (19), as:
$$\frac{{P}^{\u2033}}{P}=\frac{{\mathrm{Z\text{'}}}_{v}^{}}{{Z}_{v}}\frac{\mathrm{P\text{'}}}{P}+{Z}_{v}{Y}_{t},$$(31)where the prime symbol (′) indicates the derivative with respect to the abscissa x. This SturmLiouville equation can be rewritten into a canonical form (without the first derivative of P) by using a change of function. P is defined such as,
$$P=\psi \sqrt{{Z}_{v}},$$(32)which leads, for its first derivative, to:
$$\frac{\mathrm{P\text{'}}}{P}=\frac{\psi \text{'}}{\psi}+\frac{1}{2}\frac{{\mathrm{Z\text{'}}}_{v}^{}}{{Z}_{v}}.$$(33)
Using the second derivative of equation (32), it is obtained:
$$\frac{\mathrm{P\u2033}}{P}\frac{{\mathrm{P\text{'}}}^{2}}{{P}^{2}}=\frac{{\psi}^{\u2033}}{\psi}\frac{{\psi \mathrm{\prime}}^{2}}{{\psi}^{2}}+\frac{1}{2}\frac{{Z}_{v}^{\u2033}}{{Z}_{v}}\frac{1}{2}{\left(\frac{{\mathrm{Z\text{'}}}_{v}^{}}{{Z}_{v}}\right)}^{2}.$$(34)
Together with equation (33), this leads to:
$${\psi}^{\u2033}+{K}^{2}\left(x\right)\psi =0,$$(35)where
$${K}^{2}={Z}_{v}{Y}_{t}+\eta ,$$(36)
$$\eta =\frac{1}{2}\frac{{Z}_{v}^{\u2033}}{{Z}_{v}}\frac{3}{4}{\left(\frac{{Z\mathrm{\prime}}_{v}^{}}{{Z}_{v}}\right)}^{2}.$$(37)
Writing,
$$\frac{{\mathrm{Z\text{'}}}_{v}^{}}{{Z}_{v}}=2\frac{{\mathrm{R\text{'}}}^{}}{R}+\frac{{\mathrm{z\text{'}}}_{v}^{}}{{z}_{v}},$$(38)the coefficient η is expressed as follows:
$$\eta =\frac{{R}^{\u2033}}{R}\frac{1}{2}\frac{{z}_{v}^{\u2033}}{{z}_{v}}\frac{3}{4}\frac{{z\mathrm{\prime}}_{v}^{2}}{{z}_{v}^{2}}+\frac{\mathrm{R\text{'}}}{R}\frac{{z\mathrm{\prime}}_{v}^{}}{{z}_{v}}.$$(39)
For a cone, R″ = 0, R = x tan $\vartheta $, and R′/R = 1/x. Thus:
$$\eta =\frac{1}{2}\frac{{z}_{v}^{\u2033}}{{z}_{v}}\frac{3}{4}\frac{{z\mathrm{\prime}}_{v}^{2}}{{z}_{v}^{2}}+\frac{1}{x}\frac{{\mathrm{z\text{'}}}_{v}^{}}{{z}_{v}},$$(40)
$${z\mathrm{\prime}}_{v}^{}=2\sqrt{j}\frac{\alpha}{{x}^{2}}+6j\frac{{\alpha}^{2}}{{x}^{3}}.$$(41)
As a result, for a cone,
$$\eta =\frac{{\alpha}^{2}}{{x}^{4}}(3j+3j)=0.$$(42)
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 with high Stokes number.
3.3 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:
$${k}^{2}{\psi}^{\u2033}+q(x{)}^{2}\psi =0,$$(43)where $\psi $ is related to P by equation (32) and q is given by equation (27). This original equation is extremely simple. The term k^{−2} is the analogous to the square of the Planck constant ℏ in the literature concerning the WKB solution of the Schrödinger equation.
4 Computation of a reference solution
4.1 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}. Equation (20) give the transfer matrix of a cylindrical tube of length $\mathcal{l}$.
However the convergence with the number of frustums is more rapid if conical frustums with uniform 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, reference [6] suggests, for the frustum between x_{1} and x_{2}, to choose this radius as follows:
$$\frac{\mathcal{l}}{{R}^{\mathrm{eq}}}=\frac{{x}_{1}}{{R}_{1}}\mathrm{log}\left(1+\frac{\mathcal{l}}{{x}_{1}}\right),$$(44)where $\mathcal{l}$ = x_{2} − x_{1} is the length of the considered frustum. The origin of this expression is detailed in reference [6]. For sufficiently small conical frustums, any choice of equivalent radius R^{eq} within R_{1} and R_{2} will converge to the same solution (for small conical frustum: R^{eq} ≈ R_{1} ≈ R_{2}). Losses are computed by using equations (26) and (27) with the equivalent radius R^{eq}. The transfer matrix of a conical frustum is given by equations (15)–(17), K and Z_{c} being calculated for the the equivalent radius:
$$A=\frac{{R}_{2}}{{R}_{1}}\mathrm{cos}\left(K\mathcal{l}\right)\frac{\mathrm{sin}\left(K\mathcal{l}\right)}{K\mathcal{l}}\frac{{R}_{2}{R}_{1}}{{R}_{1}};$$(45)
$$B={Z}_{c}j\mathrm{sin}\left(K\mathcal{l}\right);$$(46)
$$D=\frac{{R}_{1}}{{R}_{2}}\mathrm{cos}\left(K\mathcal{l}\right)+\frac{\mathrm{sin}\left(K\mathcal{l}\right)}{K\mathcal{l}}\frac{{R}_{2}{R}_{1}}{{R}_{2}};$$(47)
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 St = 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 (nonradiating 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 $\mathrm{dif}\left(\mathrm{cents}\right)=1200\times {\mathrm{log}}_{2}({f}_{2}/{f}_{1})$.
Figure 2 Reference input impedance of the simplified bassoon computed by a division in 10^{5} conical frustums with ZKlosses (Bessel) computed with equations (23) and (24) for the equivalent radius ${R}_{n}^{\mathrm{eq}}$ of each conical frustum of equation (44). 
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):
$$\frac{\left\rightZ{Z}_{\mathrm{ref}}\left\right}{\left\right{Z}_{\mathrm{ref}}\left\right}=\frac{\sum _{f=20}^{f=1500}\mathrm{}\leftZ\right(f){Z}_{\mathrm{ref}}(f\left)\right}{\sum _{f=20}^{f=1500}\mathrm{}\left{Z}_{\mathrm{ref}}\right(f\left)\right}.$$(49)
Figure 3 Evolution of the norm of the relative difference with the reference impedance ${Z}_{\mathrm{ref}}$ (10^{5} conical frustums with Bessel functions) versus the length of each frustum for different approximations (cylinders with Bessel functions, conical frustums with Bessel functions, 1st or 2nd order approximation of losses). 
This error norm takes into account not only the resonance frequencies, but the entire frequency range which can influence the behaviour of the instrument. It is a usual way to observe the convergence of numerical model. The norm converges for slices of 0.1 mm 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, Figure 3 shows that the model with cylindrical frustums seems to converge also toward the reference impedance, but the convergence is slow and a too large number of cylinders would be necessary to converge (at least 10^{10}). Even if the convergence is not reached, the norm of the difference between these theoretical results and the reference impedance for the finest slicing is very small (10^{−4}).
In a musical context, the difference between the resonance peak characteristics is a more significant observation, than the error norm. The resonance frequencies f^{(i)} and magnitudes a^{(i)} of the impedance peaks are estimated by applying a second order polynomial fit on the modulus in decibel over the three samples around the maximum peaks to be more precise than the frequency step. They are computed with the finest cylindrical slicing and compared to those of the reference impedance in Figure 4. The resonance frequency deviations between the finest cylindrical slicing and the finest conical frustums, 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 three cents) (Fig. 4a). The peaks magnitude is very well estimated with a difference within 1 × 10^{−3} dB (Fig. 4b).
Figure 4 Deviation of the resonance parameters to those of the reference impedance ${Z}_{\mathrm{ref}}$ (conical frustums with ZKlosses for equivalent radius) for different approximations (cylinders with ZKlosses, conical frustums with 1st or 2nd order approximation of losses): (a) resonance frequencies and (b) resonance magnitudes. Finest slicing (10^{5} conical frustums of $2.43\times 1{0}^{5}$ m). 
4.2 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 inverse of the Stokes number (Sect. 3.1). 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. Figures 3 and 4 show that the resonance frequency deviations are within 0.1 cents for the two order of approximations of the Bessel functions (Fig. 4a). 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. 4b).
As a conclusion, the asymptotic expression of the viscothermal losses at the second order of the inverse of the Stokes number is sufficient to achieve a satisfactory accuracy.
5 WKB solution of the second order equation
We present the classical derivation of the WKB method (see Ref. [11]). The solution $\psi $ of equation (43) is sought in the following form, including an indefinite integral with respect to x:
$$\psi =g\cdot \mathrm{exp}\left(\mathrm{jk}\int u\mathrm{d}x\right),$$(50)where $g$ and u are two unknown functions and u is dimensionless. Its derivative is equal to:
$$\psi \mathrm{\prime}=g\mathrm{\prime}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+\mathrm{gjku}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right).$$(51)
Calculating the second derivative ${\psi}^{\mathrm{\u2033}}$, and inserting into equation (43), the following equation is obtained:
$${g}^{\u2033}+g{\left(\mathrm{jku}\right)}^{2}+g{k}^{2}{q}^{2}+\mathrm{jk}\left[2g\mathrm{\prime}u+\mathrm{gu}\mathrm{\prime}\right]=0.$$(52)
Because there are two unknowns functions, the vanishing of the bracket can be chosen,
$$2g\mathrm{\prime}u+\mathrm{gu}\mathrm{\prime}=0.$$(53)
This leads to the following relationship:
$$g=\frac{\mathrm{\Lambda}}{\sqrt{u}},$$(54)where Λ is a constant. The other part of equation (52) yields:
$$\frac{{g}^{\u2033}}{g}={k}^{2}\left({u}^{2}{q}^{2}\right).$$(55)
If u is a solution, −u is also a solution. Therefore $\psi $ is the superposition of two solutions, as follows:
$$\psi ={u}^{1/2}\left[{\psi}^{+}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+{\psi}^{}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)\right],$$(56)where ${\psi}^{\pm}$ are two constants. u is the unique unknown function. It is the solution of the following equation, which is derived from equations (54) and (55):
$${k}^{2}{u}^{2}({u}^{2}{q}^{2})=\frac{3}{4}{\mathrm{u\text{'}}}^{2}\frac{1}{2}{u}^{\u2033}u.$$(57)
For infinite k, a trivial solution is u = q. For large k, the solution is sought in the form of an asymptotic expansion:
$$u=q+v/{k}^{2}+w/{k}^{4}+\dots $$(58)
Simplifying by the factor 1/k^{2}, the term v can be obtained:
$$\left(v+\frac{w}{{k}^{2}}\right)\left(u+q\right){u}^{2}=\frac{3}{4}{u\text{'}}^{2}\frac{1}{2}{u}^{\u2033}u.$$(59)therefore, at the second order of 1/k :
$$v=3\frac{{\mathrm{q\text{'}}}^{2}}{8{q}^{3}}\frac{{q}^{\u2033}}{4{q}^{2}}.$$(60)
Another useful expression for the function v is:
$$v=\frac{1}{4}{\left[\frac{\mathrm{q\text{'}}}{{q}^{2}}\right]}^{\text{'}}\frac{{q}^{\mathrm{\prime}2}}{8{q}^{3}}.$$(61)
The expression of the fourth order coefficient leads to complicated results:
$$w=\frac{5}{2}\frac{{v}^{2}}{q}+\frac{3}{4}\frac{\mathrm{q\text{'}v\text{'}}}{{q}^{3}}\frac{1}{4}\left[\frac{\mathrm{q\u2033}v}{{q}^{3}}+\frac{\mathrm{v\text{'}\text{'}}}{{q}^{2}}\right].$$(62)This expression can be a starting point for further studies. However, the challenge of the present paper being to avoid large complication in the expressions of the transfer matrix, this expression is not used throughout the rest of this study.
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 (Sect. 2):
$$\widehat{P}={\psi}^{+}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+{\psi}^{}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right).$$(63)
Using equations (32) and (56), the relationship between the pressure P and $\widehat{P}$ is found to be:
$$P=\widehat{P}\sqrt{\frac{{Z}_{v}}{u}},$$(64)thus,
$$\frac{\mathrm{P\text{'}}}{P}=\frac{\widehat{P}\text{'}}{\widehat{P}}+\frac{1}{2}\left[\frac{{\mathrm{z\text{'}}}_{v}^{}}{{z}_{v}}\frac{\mathrm{u\text{'}}}{u}\right].$$(65)
The variable $\widehat{U}$ is defined as:
$$\widehat{U}={\psi}^{+}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right){\psi}^{}\mathrm{exp}\left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right).$$(66)Therefore,
$$\widehat{U}=j\widehat{P}\mathrm{\prime}/\left(\mathrm{ku}\right).$$(67)
For the variables $\widehat{P}$ and $\widehat{U}$, thanks to equations (63) and (66) the following transfer matrix relationship can be introduced:
$${\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{1}=\left(\begin{array}{ll}\mathrm{cos}X& j\mathrm{sin}X\\ j\mathrm{sin}X& \mathrm{cos}X\end{array}\right){\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{2},$$(68)with
$$X=k{\int}_{{x}_{1}}^{{x}_{2}}\mathrm{}u\mathrm{d}x.$$(69)
The development of the expression of u makes appearing the integration of 1/R (Eqs. (29) and (58)), which is coherent with the equivalent radius R^{eq} chosen in equation (44). The relationship between (P, U) and ($\widehat{P}$, $\widehat{U}$) comes from equations (18) and (64):
$$U=\frac{\mathrm{P\text{'}}}{{Z}_{v}}=\frac{\mathrm{P\text{'}}}{P}\frac{\widehat{P}}{\sqrt{{Z}_{v}u}}.$$(70)Thus, with equations (65) and (67):
$$U=\frac{1}{\sqrt{{Z}_{v}u}}\left[\mathrm{jku}\widehat{U}\frac{1}{2}\left(\frac{{\mathrm{Z\text{'}}}_{v}^{}}{{Z}_{v}}\frac{\mathrm{u\text{'}}}{u}\right)\widehat{P}\right].$$(71)
The transition matrix at each extremity of the truncated cone is derived from equations (64) and (71) as follows:
$${\left(\begin{array}{l}P\\ U\end{array}\right)}_{i}=\left(\begin{array}{ll}{a}_{i}& 0\\ {c}_{i}& {d}_{i}\end{array}\right){\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{i},$$(72)with,
$$\begin{array}{c}{a}_{i}=\sqrt{{Z}_{\mathrm{vi}}{u}_{i}}\\ {d}_{i}=\mathrm{jk}\sqrt{{u}_{i}{Z}_{\mathrm{vi}}}\\ {c}_{i}=\frac{1}{2}\frac{1}{\sqrt{{Z}_{\mathrm{vi}}{u}_{i}}}\left(\frac{{\mathrm{Z\text{'}}}_{{v}_{i}}^{}}{{Z}_{{v}_{i}}}\frac{{\mathrm{u\text{'}}}_{i}^{}}{{u}_{i}}\right).\end{array}$$
The determinant of the matrix is ${a}_{i}{d}_{i}=\mathrm{jk}$ and is independent of the extremity (1 or 2). The inverse matrix is:
$${\left(\begin{array}{l}\widehat{P}\\ \widehat{U}\end{array}\right)}_{i}=\left(\begin{array}{ll}1/{a}_{i}& 0\\ {c}_{i}/\left({a}_{i}{d}_{i}\right)& 1/{d}_{i}\end{array}\right){\left(\begin{array}{l}P\\ U\end{array}\right)}_{i}.$$(73)
The product of the three matrices between abscissae 1 and 2 gives the final transfer matrix:
$${\left(\begin{array}{l}P\\ U\end{array}\right)}_{1}=\left(\begin{array}{ll}A& B\\ C& D\end{array}\right){\left(\begin{array}{l}P\\ U\end{array}\right)}_{2};$$(74)with,
$$\begin{array}{ll}A& =\frac{{a}_{1}}{{a}_{2}}\left[\mathrm{cos}X\frac{{c}_{2}}{{d}_{2}}j\mathrm{sin}X\right]\\ D& =\frac{{d}_{1}}{{d}_{2}}\left[\mathrm{cos}X+\frac{{c}_{1}}{{d}_{1}}j\mathrm{sin}X\right]\\ B& =\frac{{a}_{1}}{{d}_{2}}j\mathrm{sin}\left(X\right).\\ & \end{array}$$(75)
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:
$$A=\sqrt{\frac{{Z}_{v1}}{{Z}_{v2}}}\sqrt{\frac{{u}_{2}}{{u}_{1}}}\left[\mathrm{cos}X+\frac{1}{2}\frac{1}{k{u}_{2}}\left(\frac{{\mathrm{Z\text{'}}}_{v2}^{}}{{Z}_{v2}}\frac{{\mathrm{u\text{'}}}_{2}^{}}{{u}_{2}}\right)\mathrm{sin}X\right];$$(76)
$$D=\sqrt{\frac{{Z}_{v2}}{{Z}_{v1}}}\sqrt{\frac{{u}_{1}}{{u}_{2}}}\left[\mathrm{cos}X\frac{1}{2}\frac{1}{k{u}_{1}}\left(\frac{{\mathrm{Z\text{'}}}_{v1}^{}}{{Z}_{v1}}\frac{{\mathrm{u\text{'}}}_{1}^{}}{{u}_{1}}\right)\mathrm{sin}X\right];$$(77)
$$B=\frac{1}{k}\mathrm{sin}X\sqrt{\frac{{Z}_{v1}{Z}_{v2}}{{u}_{1}{u}_{2}}};$$(78)
$$C=\frac{\mathrm{AD}1}{B}.$$(79)
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 equation (58), we need the equivalent expressions for the functions q and v, their integral and their derivative.
6.2 Expansions of q and v with respect to the Stokes number
The starting point is the expansion of the quantity q, as written in equation (27), at the second order of the inverse of the Stokes number. The indefinite integral ${I}_{q}=k\int \mathrm{}q\mathrm{d}x$ and the derivative q′ are given by:
$${I}_{q}=x+\gamma \alpha \mathrm{ln}\left(x\right)\delta \frac{{\alpha}^{2}}{x};$$(80)
$$q\text{'}=\gamma \frac{\alpha}{{x}^{2}}2\delta \frac{{\alpha}^{2}}{{x}^{3}}.$$(81)
Notice that x is a variable with dimension, thus the natural 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 equation (61). The integral of the second term of the expression (61) is complicated, but a numerical evaluation is found to be very small, therefore,
$${I}_{v}\simeq \frac{1}{4}\frac{q\text{'}}{{q}^{2}}.$$(82)
In order to calculate v′, we need the calculation of the second and third derivatives of q′ (Eq. (60):
$$\begin{array}{c}{q}^{\u2033}=\frac{2\gamma \alpha}{{x}^{3}}+\frac{6\delta {\alpha}^{2}}{{x}^{4}}\\ {q}^{\mathrm{\u2034}}=\frac{6\gamma \alpha}{{x}^{4}}+\frac{24\delta {\alpha}^{2}}{{x}^{5}}\\ v\text{'}=\frac{1}{{q}^{3}}\left[\frac{5}{4}q\mathrm{\text{'}q\text{'}\text{'}}\frac{9}{8}\frac{{q}^{\text{'}3}}{q}\frac{q{q}^{\mathrm{\u2034}}}{4}\right].\end{array}$$(83)
Finally, we need to successively compute the following expressions:
Compute α, γ, δ with equations (30) and (28) then ${z}_{v},{z\text{'}}_{v}^{},q,{I}_{q},q\mathrm{\prime},\mathrm{q\u2033}$ and $\mathrm{q\u2034}$ with equations (29), (41), (27), (80), (81), and (83).
Compute Z_{v} and ${Z\text{'}}_{v}^{}/{Z}_{v}$ with equations (22) and (38).
Compute v′ and I_{v} then v with equations (82), (83), and (60).
Compute u and u′ with equation (58).
7 Comparison of the zeroth and second orders of the WKB solution
7.1 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 Section 4.1.
We first notice that the obtained transfer matrix given by equations (76)–(78) 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.
Because the solution from the WKB method is aimed at limiting the number of frustrums to only one, using a slicing is not consistent with this aim, and therefore not useful. This has to be distinguished from the slicing with an approximate equation (based upon averaged losses). More explicitly, when no losses are present, the slicing of a cone does not improve the result given by one matrix, because the equation to be solved is exact. Comparing the classical method of slicing in cylindrical or conical segment (with averaged losses), it can be noticed that the classical method converges when the number of frustums tends to infinity, while the WKB method could converge by extending the order of expansion of the function u. This convergence is out of the scope of the paper.
The analytical limit of validity of the WKB solution can be given by the comparison between the terms due to the viscothermal effects and their variations with the radius. At the first order of the inverse of the Stokes number and at the second order of the WKB solution, the function u is written as:
$$u=1+\frac{\gamma \alpha}{x}\left[1\frac{1}{2{k}^{2}{x}^{2}}\right].$$(84)
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:
$${I}_{u}=x+\gamma \alpha \left[\mathrm{ln}\left(x\right)\frac{1}{4{k}^{2}{x}^{2}}\right],$$(85)and,
$$u\text{'}=\frac{\gamma \alpha}{{x}^{2}}\left[1\frac{3}{2{k}^{2}{x}^{2}}\right].$$(86)
Using equation (85) the quantity X (Eq. (69)) can be calculated. The result is:
$$X=k\mathcal{l}+\frac{\gamma}{\mathrm{tan}\vartheta}\sqrt{\frac{c{\mathcal{l}}_{v}}{\omega}}\left[\mathrm{ln}\left\{1+\frac{\mathcal{l}}{{x}_{1}}\right\}\frac{1}{4{k}^{2}{x}_{2}^{2}}+\frac{1}{4{k}^{2}{x}_{1}^{2}}\right].$$(87)
Therefore, an approximate condition of validity can be estimated: x_{1} > 1. For the first peak of a divergent cone, kx_{1} ≈ π/x_{2}, thus the condition can be written as:
$$\pi \frac{{x}_{1}}{{x}_{2}}=k{x}_{1}>1,$$or $\mathcal{l}<(\pi 1){x}_{1}$. This is not satisfied for the first peak of a bassoon. However, above the third or fourth peak, this becomes acceptable.
7.2 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. 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 (Fig. 5). Both correspond to frequencies for which kx_{1} < 1. However, the zeroth order approximation seems better at these peaks (Fig. 5). This observation can be quantified by applying the error norm of equation (49) 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).
Figure 5 Input impedance of the simplified bassoon (reference impedance) compared with the WKB formulation on one conical frustums at the second and zeroth order (in both WKB solutions, the second order order is used for the Stokes number). 
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. 6a) the zeroth order is better. It induces a deviation about −30 cents and −4 dB 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. 6b), the second order is slightly better, even if both orders show acceptable deviations to the reference impedance (<1 cents and <0.5 dB).
Figure 6 Difference between the resonance parameters of the WKB simulations (second and zeroth order) to those of the reference impedance: top: frequencies and bottom: magnitudes. (a) Low frequencies and (b) high frequencies (the vertical scales are different). In both WKB solutions, the second order approximation is used for the viscothermal effects. 
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.
8 A simplified, general formula
8.1 Plane waves
It can be useful to search for a further simplification of the zeroth order solution. We start from equations (76)–(78) with q_{1} = u_{1} and q_{2} = u_{2}. First of all, we notice the following relationship: ${Z}_{v}/q=\mathrm{jk}{Z}_{c}$. The transfer matrix can be written as follows:
$$A=\sqrt{\frac{{Z}_{c1}}{{Z}_{c2}}}\left[\mathrm{cos}X+\frac{1}{2}\frac{1}{k{q}_{2}}\frac{{\mathrm{Z\text{'}}}_{c2}^{}}{{Z}_{c2}}\mathrm{sin}X\right];$$(88)
$$A=\sqrt{\frac{{Z}_{c2}}{{Z}_{c1}}}\left[\mathrm{cos}X\frac{1}{2}\frac{1}{k{q}_{1}}\frac{{\mathrm{Z\text{'}}}_{c1}^{}}{{Z}_{c1}}\mathrm{sin}X\right];$$(89)
$$B=\frac{1}{k}\mathrm{sin}X\sqrt{\frac{{Z}_{v1}{Z}_{v2}}{{u}_{1}{u}_{2}}}.$$(90)
At the second order of the inverse of the Stokes number, the characteristic impedance is written as:
$${Z}_{c}=\frac{\rho c}{S}\left[1+0.37\sqrt{2j}S{t}^{1}1.147\mathrm{jS}{t}^{2}\right].$$(91)
Therefore, by comparing with the expression of the wavenumber K = kq (see Eq. (29)), the viscothermal effects intervene significantly more in K than in $\sqrt{{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:
$$A=\frac{{R}_{2}}{{R}_{1}}\mathrm{cos}X\frac{1}{k{q}_{2}{x}_{1}}\mathrm{sin}X;$$(92)
$$D=\frac{{R}_{1}}{{R}_{2}}\mathrm{cos}X+\frac{1}{k{q}_{1}{x}_{2}}\mathrm{sin}X;$$(93)
$$B=\frac{\rho c}{\pi {R}_{1}{R}_{2}}j\mathrm{sin}X.$$(94)
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 viscothermal effects. 
This approximation differs from other approximations published by some authors. Concerning the original approach by Kulik [12], there are two differences with the present approach: firstly the series impedance Z_{v} does not vary with the radius (${q}_{1}={q}_{2}=1$ in our Eqs. (89) and (90)); secondly the differential equation to be solved is not explicitly written: the starting point built with the equations (2) and (8) is not clearly justified, even if the quantity $\overline{k}L$ is the same as our quantity X. Other approaches were given for the calculation of the impedance peaks [6, 13] or of the resonance frequencies [14], the results being also consistent with equation (69), which is based upon the calculation of a mean value of the inverse of the radius. The approach by Nederveen [14] was a perturbation method, limited to the first order of the variation of losses with the radius, and the computation was limited to the resonance frequencies (real wavenumbers). The generalization to complex wavenumbers could be probably possible.
We emphasize that this paper considers an extreme case, the bassoon. Figure 8 shows that for a wider instrument, such as a simplified baritone saxophone (R_{1} = 6.75 mm, R_{2} = 60 mm, $\mathcal{l}$ = 2.38 m), the final results are much better for the first peaks. Furthermore the convergence study to the reference result justifies a classical approximation of viscothermal effects in cones: the second order of the asymptotic expansion with respect to the inverse of the Stokes number does not significantly improve the accuracy of the computation. This was not clear in [13].
Figure 8 Difference between the resonance parameters of the WKB simulations at the zeroth order (complete and simplified) to those of a reference impedance for a simplified baritone saxophone: top: frequencies and bottom: magnitudes. In both simulations, the second order approximation is used for the viscothermal effects. 
8.2 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 [6]. For such cones, viscothermal effects are weak, and the WKB solution of zeroth order is satisfactory. It is not necessary to repeat the complete analysis. We propose the following formula, which is a modification of equations (8)–(10):
$$\begin{array}{c}A=\frac{{r}_{2}}{{r}_{1}}\mathrm{cos}{X}_{s}\frac{\mathrm{sin}{X}_{s}}{k{q}_{2}{r}_{1}};\\ D=\frac{{r}_{1}}{{r}_{2}}\mathrm{cos}{X}_{s}+\frac{\mathrm{sin}{X}_{s}}{k{q}_{1}{r}_{2}};\\ \begin{array}{c}B=\frac{\rho c}{{\mathrm{\Sigma}}_{2}}\frac{{r}_{2}}{{r}_{1}}j\mathrm{sin}{X}_{s};\\ C=\frac{\mathrm{AD}1}{B};\end{array}\end{array}$$(95)with,
$${X}_{s}=k\left({I}_{q2}{I}_{q1}\right),$$(96)and,
$${I}_{q}=r+\gamma \alpha \mathrm{ln}\left(r\right)\delta \frac{{\alpha}^{2}}{r}.$$(97)
9 Conclusion
To our knowledge the propagation equation in cones with viscothermal effects was not derived in previous papers. In the frequency domain, at the second order of the asymptotic expansion with respect to the inverse of 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 will 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 viscothermal 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 (95) are a satisfactory compromise for conical instruments.
The aim of the present work is the analytical derivation of the WKB solution, in order to avoid the slicing in frustums. This derivation is done at the second order of the asymptotic expansion with respect the inverse of the Stokes number, and is shown to be sufficiently accurate for the application to low wind instruments. For other applications, to extend the expansion to other orders, or even to consider the very small Stokes number (capillary tubes) could be done by using the other asymptotic expansion (see [9, 15]).
For real woodwind instruments, tone holes or change of conicity impose to slice the main bore in several conical frustums. For the narrow parts where the condition $k{x}_{1}>1$ is not respected, it is better to use the usual method (slice in several frustums with equivalent losses), but for the wide parts, the WKB solution can improve and ease the computation of the transfer matrices.
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 [6], but the corresponding work should be heavy.
Acknowledgments
This work has been partly supported by the french Agence Nationale de la Recherche (ANR16LCV2000701 Liamfi project), in cooperation with BuffetCrampon. We thank JeanPierre Dalmont for useful discussions.
Conflict of interest
Author declared no conflict of interests.
References
 G.R. Plitnik, W.J. Strong: Numerical method for calculating input impedances of the oboe. Journal of the Acoustical Society of America 65 (1979) 816–825. [CrossRef] [Google Scholar]
 D.H. Keefe: Woodwind tonehole acoustics and the spectrum transformation function. Ph.D. Thesis. CaseWestern Reserve University, Cleveland, OH, 1981. [Google Scholar]
 R. Caussé, J. Kergomard, X. Lurton: Input impedance of brass musical instruments. Journal of the Acoustical Society of America 75 (1984) 241–254. [CrossRef] [Google Scholar]
 P. Eveno, J.P. Dalmont, R. Caussé, J. Gilbert: Wave propagation and radiation in a horn: Comparisons between models and measurements. Acta Acustica United With Acustica 98 (2012) 158–165. [CrossRef] [Google Scholar]
 C. Zwikker, C.W. Kosten: Sound Absorbing Materials. Elsevier, New York, 1949. [Google Scholar]
 A. Chaigne, J. Kergomard: Acoustics of Musical Instruments. New York, USA, 2016. [CrossRef] [Google Scholar]
 N. Fletcher, T. Rossing: Physics of Musical Instruments. 2nd ed., Springer, New York, USA, 1998. [Google Scholar]
 D.H. Keefe: Acoustical wave propagation in cylindrical ducts: Transmission line parameter approximations for isothermal and nonisothermal boundary conditions. Journal of the Acoustical Society of America 75 (1984) 58–62. [CrossRef] [Google Scholar]
 T. Hélie: Unidimensional models of acoustic propagation in axisymmetric waveguides. Journal of the Acoustical Society of America 114 (2003) 2633–2647. [CrossRef] [PubMed] [Google Scholar]
 T. Hélie, T. Hézard, R. Mignot, D. Matignon: Onedimensional acoustic models of horns and comparison with measurements. Acta Acustica United With Acustica 99 (2013) 960–974. [CrossRef] [Google Scholar]
 A. Voros: The return of the quartic oscillator. The complex WKB method. Annales de l’I.H.P. Physique Théorique 39 (1983) 211–338. [Google Scholar]
 Y. Kulik: Transfer matrix of conical waveguides with any geometric parameters for increased precision in computer modeling. Journal of the Acoustical Society of America 122 (2007) EL 179–EL 184. [CrossRef] [Google Scholar]
 J. Kergomard: Quasistanding waves in horns with viscothermal losses at the walls: computation of the impedance. Acustica 48 (1981) 31–43. (in French). [Google Scholar]
 C.J. Nederveen: Acoustical Aspects Woodwind Instruments. 2nd ed., Northern Illinois University Press, Dekalb, IL, 1998. [Google Scholar]
 J. Kergomard, R. Caussé: Measurement of acoustic impedance using a capillary: An attempt to achieve optimization. Journal of the Acoustical Society of America 79 (1986) 1129–1140. [CrossRef] [Google Scholar]
Cite this article as: Ernoult A & Kergomard J. 2020. Transfer matrix of a truncated cone with viscothermal losses: application of the WKB method. Acta Acustica, 4, 7.
All Figures
Figure 1 Sketch of a conical tube. $L={r}_{2}{r}_{1}=\mathcal{l}/\mathrm{cos}\vartheta .$ 

In the text 
Figure 2 Reference input impedance of the simplified bassoon computed by a division in 10^{5} conical frustums with ZKlosses (Bessel) computed with equations (23) and (24) for the equivalent radius ${R}_{n}^{\mathrm{eq}}$ of each conical frustum of equation (44). 

In the text 
Figure 3 Evolution of the norm of the relative difference with the reference impedance ${Z}_{\mathrm{ref}}$ (10^{5} conical frustums with Bessel functions) versus the length of each frustum for different approximations (cylinders with Bessel functions, conical frustums with Bessel functions, 1st or 2nd order approximation of losses). 

In the text 
Figure 4 Deviation of the resonance parameters to those of the reference impedance ${Z}_{\mathrm{ref}}$ (conical frustums with ZKlosses for equivalent radius) for different approximations (cylinders with ZKlosses, conical frustums with 1st or 2nd order approximation of losses): (a) resonance frequencies and (b) resonance magnitudes. Finest slicing (10^{5} conical frustums of $2.43\times 1{0}^{5}$ m). 

In the text 
Figure 5 Input impedance of the simplified bassoon (reference impedance) compared with the WKB formulation on one conical frustums at the second and zeroth order (in both WKB solutions, the second order order is used for the Stokes number). 

In the text 
Figure 6 Difference between the resonance parameters of the WKB simulations (second and zeroth order) to those of the reference impedance: top: frequencies and bottom: magnitudes. (a) Low frequencies and (b) high frequencies (the vertical scales are different). In both WKB solutions, the second order approximation is used for the viscothermal effects. 

In the text 
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 viscothermal effects. 

In the text 
Figure 8 Difference between the resonance parameters of the WKB simulations at the zeroth order (complete and simplified) to those of a reference impedance for a simplified baritone saxophone: top: frequencies and bottom: magnitudes. In both simulations, the second order approximation is used for the viscothermal effects. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.