© A. Ernoult and J. Kergomard, Published by EDP Sciences, 2020

1 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 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 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 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₁, output radius R₂ (Fig. 1). 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₂ − r₁ is equal to (R₂ − R₁)/sin $\vartheta$. The area of the spherical wavefront is $\mathrm{\Sigma }=2\pi r^2\left[1-\cos \vartheta \right]=2\pi R^2/\left[1+\cos \vartheta \right]$.

Figure 1 Sketch of a conical tube. $L=r_2-r_1=\mathcal{l}/\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₂ − x₁ is equal to L cos $\vartheta$. The area of the planar wavefront is simply S = πR². 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}\cos \mathrm{kL}& j\sin \mathrm{kL}\\ j\sin \mathrm{kL}& \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}\cos \left(\mathrm{kL}\right)-\frac{\sin \left(\mathrm{kL}\right)}{k{r}_1};$$(8)

$$D=\frac{r_1}{r_2}\cos \left(\mathrm{kL}\right)+\frac{\sin \left(\mathrm{kL}\right)}{k{r}_2};$$(9)

$$B=\frac{\varrho c}{{\mathrm{\Sigma }}_2}\frac{r_2}{r_1}j\sin \left(\mathrm{kL}\right);$$(10)

$$C=(\mathrm{AD}-1)/B.$$(11)

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}\cos \left(\mathrm{kL}\right)-\frac{\sin \left(\mathrm{kL}\right)}{\mathrm{kL}}\frac{R_2-R_1}{R_1};$$(12)

$$D=\frac{R_1}{R_2}\cos \left(\mathrm{kL}\right)+\frac{\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+\cos \vartheta }{2}j\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}\cos \left(k\mathcal{l}\right)-\frac{\sin \left(k\mathcal{l}\right)}{k\mathcal{l}}\frac{R_2-R_1}{R_1};$$(15)

$$D=\frac{R_1}{R_2}\cos \left(k\mathcal{l}\right)+\frac{\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\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=\cos K\mathcal{l}\\ B=j Z_c\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\right|R$. 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 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, 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\tan \vartheta \left|k_v\right|$,

$$\alpha =\frac{x}{\mathrm{St}}=\frac{1}{\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^{″}}{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 Sturm-Liouville 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″}}{P}-\frac{{\mathrm{P\text{'}}}^2}{P^2}=\frac{{\psi }^{″}}{\psi }-\frac{{\psi \mathrm{\prime }}^2}{{\psi }^2}+\frac{1}{2}\frac{Z_v^{″}}{Z_v}-\frac{1}{2}{\left(\frac{{\mathrm{Z\text{'}}}_v^{}}{Z_v}\right)}^2.$$(34)

Together with equation (33), this leads to:

$${\psi }^{″}+K^2\left(x\right)\psi =0,$$(35)where

$$K^2=-Z_v{Y}_t+\eta ,$$(36)

$$\eta =\frac{1}{2}\frac{Z_v^{″}}{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^{″}}{R}\frac{1}{2}\frac{z_v^{″}}{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^{″}}{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 }^{″}+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⁻² 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₁ and x₂, to choose this radius as follows:

$$\frac{\mathcal{l}}{R^{\mathrm{eq}}}=\frac{x_1}{R_1}\log \left(1+\frac{\mathcal{l}}{x_1}\right),$$(44)where $\mathcal{l}$ = x₂ − x₁ 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₁ and R₂ will converge to the same solution (for small conical frustum: R^eq ≈ R₁ ≈ R₂). 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}\cos \left(K\mathcal{l}\right)-\frac{\sin \left(K\mathcal{l}\right)}{K\mathcal{l}}\frac{R_2-R_1}{R_1};$$(45)

$$B=Z_{c}j\sin \left(K\mathcal{l}\right);$$(46)

$$D=\frac{R_1}{R_2}\cos \left(K\mathcal{l}\right)+\frac{\sin \left(K\mathcal{l}\right)}{K\mathcal{l}}\frac{R_2-R_1}{R_2};$$(47)

$$C=(\mathrm{AD}-1)/B.$$(48)

The example of a simplified bassoon (R₁ = 0.002 m; R₂ = 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 (non-radiating open pipe) (Fig. 2). The computation is carried out for the frequency range [20, 10⁴] 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 {\log }_2(f_2/f_1)$.

Figure 2 Reference input impedance of the simplified bassoon computed by a division in 105 conical frustums with ZK-losses (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|\right|Z-Z_{\mathrm{ref}}\left|\right|}{\left|\right|Z_{\mathrm{ref}}\left|\right|}=\frac{\sum _{f=20}^{f=1500}\mathrm{}\left|Z\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}}$ (105 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⁴ frustums for the considered tube. With this slicing the variation of the radius for one frustum is about 7 × 10⁻⁷ m. Throughout this paper, the reference impedance Z_ref corresponds to that computed with the finest slicing (10⁵ conical frustums of 2.43 × 10⁻⁵ 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¹⁰). 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⁻⁴).

In a musical context, the difference between the resonance peak characteristics is a more significant observation, than the error norm. The resonance frequencies f⁽ⁱ⁾ and magnitudes a⁽ⁱ⁾ 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⁻³ dB (Fig. 4b).

Figure 4 Deviation of the resonance parameters to those of the reference impedance $Z_{\mathrm{ref}}$ (conical frustums with ZK-losses for equivalent radius) for different approximations (cylinders with ZK-losses, conical frustums with 1st or 2nd order approximation of losses): (a) resonance frequencies and (b) resonance magnitudes. Finest slicing (105 conical frustums of $2.43\times 10^{-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⁻² for the first order and 3 × 10⁻⁴ 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 visco-thermal 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 \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 }\exp \left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+\mathrm{gjku} \exp \left(\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right).$$(51)

Calculating the second derivative ${\psi }^{\mathrm{″}}$, and inserting into equation (43), the following equation is obtained:

$$g^{″}+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^{″}}{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 }^{+} \exp \left(-\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+{\psi }^{-} \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}^{″}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², 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}^{″}u.$$(59)therefore, at the second order of 1/k :

$$v=3\frac{{\mathrm{q\text{'}}}^2}{8q^3}-\frac{q^{″}}{4q^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}}{8q^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″}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 }^{+} \exp \left(-\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)+{\psi }^{-} \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 }^{+} \exp \left(-\mathrm{jk}\int \mathrm{}u\mathrm{d}x\right)-{\psi }^{-}\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}\cos X& j\sin X\\ j\sin X& \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[\cos X-\frac{c_2}{d_2}j\sin X\right]\\ D& =\frac{d_1}{d_2}\left[\cos X+\frac{c_1}{d_1}j\sin X\right]\\ B& =\frac{a_1}{d_2}j\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[\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)\sin X\right];$$(76)

$$D=\sqrt{\frac{Z_{v2}}{Z_{v1}}}\sqrt{\frac{u_1}{u_2}}\left[\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)\sin X\right];$$(77)

$$B=\frac{1}{k}\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 \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₂/x₁) 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}^{″}=\frac{2\gamma \alpha }{x^3}+\frac{6\delta {\alpha }^2}{x^4}\\ q^{\mathrm{‴}}=-\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{‴}}}{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″}$ and $\mathrm{q‴}$ 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).

• Compute I_u and X with equations (58) and (69).

• Finally compute A, B, C, D with equations (76)–(79).

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⁻¹. 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 visco-thermal 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}{2k^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[\ln \left(x\right)-\frac{1}{4k^2{x}^2}\right],$$(85)and,

$$u\text{'}=-\frac{\gamma \alpha }{x^2}\left[1-\frac{3}{2k^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 }{\tan \vartheta }\sqrt{\frac{c{\mathcal{l}}_v}{\omega }}\left[\ln \left\{1+\frac{\mathcal{l}}{x_1}\right\}-\frac{1}{4k^2{x}_2^2}+\frac{1}{4k^2{x}_1^2}\right].$$(87)

Therefore, an approximate condition of validity can be estimated: x₁ > 1. For the first peak of a divergent cone, kx₁ ≈ π/x₂, 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⁵ 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. 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₁ > 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₁ in Figure 6. At low frequencies (kx₁ < 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₁ > 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 visco-thermal 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₁ = u₁ and q₂ = u₂. 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[\cos X+\frac{1}{2}\frac{1}{k{q}_2}\frac{{\mathrm{Z\text{'}}}_{c2}^{}}{Z_{c2}}\sin X\right];$$(88)

$$A=\sqrt{\frac{Z_{c2}}{Z_{c1}}}\left[\cos X-\frac{1}{2}\frac{1}{k{q}_1}\frac{{\mathrm{Z\text{'}}}_{c1}^{}}{Z_{c1}}\sin X\right];$$(89)

$$B=\frac{1}{k}\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 visco-thermal 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₁ and x₂. As a consequence, the following simplified formula can be useful:

$$A=\frac{R_2}{R_1}\cos X-\frac{1}{k{q}_2{x}_1}\sin X;$$(92)

$$D=\frac{R_1}{R_2}\cos X+\frac{1}{k{q}_1{x}_2}\sin X;$$(93)

$$B=\frac{\rho c}{\pi R_1{R}_2}j\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 visco-thermal 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₁ = 6.75 mm, R₂ = 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 visco-thermal 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, visco-thermal 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}\cos X_s-\frac{\sin X_s}{k{q}_2{r}_1};\\ D=\frac{r_1}{r_2}\cos X_s+\frac{\sin X_s}{k{q}_1{r}_2};\\ \begin{array}{c}B=\frac{\rho c}{{\mathrm{\Sigma }}_2}\frac{r_2}{r_1}j\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 \ln \left(r\right)-\delta \frac{{\alpha }^2}{r}.$$(97)

9 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 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₁ 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 visco-thermal effects in the coefficients q₁ and q₂ (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 (ANR16-LCV2-0007-01 Liamfi project), in cooperation with Buffet-Crampon. We thank Jean-Pierre Dalmont for useful discussions.

Conflict of interest

Author declared no conflict of interests.

References

All Figures

	Figure 1 Sketch of a conical tube. $L=r_2-r_1=\mathcal{l}/\cos \vartheta .$
In the text

	Figure 2 Reference input impedance of the simplified bassoon computed by a division in 105 conical frustums with ZK-losses (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}}$ (105 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 ZK-losses for equivalent radius) for different approximations (cylinders with ZK-losses, conical frustums with 1st or 2nd order approximation of losses): (a) resonance frequencies and (b) resonance magnitudes. Finest slicing (105 conical frustums of $2.43\times 10^{-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 visco-thermal 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 visco-thermal 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