Alexis Thibault
email: alexis.thibault@inria.fr

Juliette Chabassier

Henri Boutin

Thomas Hélie

Transmission line coefficients for viscothermal acoustics in conical tubes

Keywords: viscothermal acoustics, transmission line model, conical tubes, hydraulic radius

Viscothermal acoustic propagation in gases contained in rigid straight or conical tubes is considered. Under the assumption that the wavelength is much larger than both the boundary layer thickness and the tube radius, pressure and flow are shown to be solutions of a pair of coupled 1D differential equations, formulated as transmission line equations involving complex loss coefficients. The derivation of these loss coefficients, which is usually accomplished in cylinders, is generalized here to conical geometries. In the well-kown case of circular cylinders, the Zwikker-Kosten (ZK) theory is recovered. For circular cones, the expression of the loss coefficients is derived. It involves complex-order spherical harmonics, instead of Bessel functions for circular cylinders, and makes the hydraulic radius appear as a natural relevant geometrical parameter.

We show that replacing the classical radius by the hydraulic radius in the ZK theory provides an affordable and accurate approximation of the analytic model derived for cones. The proposed formulas are used to compute the input impedance of a cone, and compared with a 3D reference.

In an ideal setting, using the spherical harmonics or the hydraulic radius in the 1D method accurately approximates the full 3D method, and allows to increase accuracy by approximately two orders of magnitude compared to the ZK theory.

Thermal conductivity κ 2.57 × 10 -2 J m -1 s -1 K -1

Specific heat w/ constant pressure Cp 1004 J kg -1 K -1

Specific gas constant R0

γ-1 γ Cp 288 J kg -1 K -1

Half opening angle of the cone Θ 0 to 0. Boundary-layer-wavelength ratio kR/ Sh µω/(c 2 0 ρ0) 10 -4 to 4 × 10 -3 (no unit)

Table 1: Notations used throughout this paper. Four characteristic lengths are exhibited. The four significant unitless coefficients from [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF] are highlighted. Values of most physical coefficients are calculated from [2, (5.142)]

for air at 20 °C, and other values are estimated from them. Bulk viscosity is estimated from that of N 2 [3, p.9].

.

Introduction

In the context of musical acoustics, virtual instrument prototyping requires accurate models of sound propagation in tubes with varying cross-section. Such models are commonly formulated as transmission line equations, a pair of coupled partial differential equations with a single spatial dimension [START_REF] Rienstra | Webster's horn equation revisited[END_REF][START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF][START_REF] Hélie | Unidimensional models of acoustic propagation in axisymmetric waveguides[END_REF]. In order to account for the presence of boundary layers near the tube walls, the coefficients of these equations should depend on the viscosity and thermal properties of the fluid. One possible expression for the transmission line coefficients of cylindrical tubes was initially proposed by Zwikker and Kosten [START_REF] Zwikker | Sound Absorbing Materials[END_REF], and is now widely used to model tubes of any shape [START_REF] Caussé | Input impedance of brass musical instruments -comparison between experiment and numerical models[END_REF][START_REF] Bilbao | Direct simulation of reed wind instruments[END_REF][START_REF] Chaigne | Acoustics of musical instruments[END_REF][START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF]. The current work aims to question whether the use of these coefficients is justified in tubes with varying cross-section in spite of the fact that they are derived for cylinders. This issue is addressed by considering the special case of conical tubes.

The motion of air in a thin tube of varying section can be approximated by the following "horn equation" [START_REF] Rienstra | Webster's horn equation revisited[END_REF], which relates the spatial evolution of the acoustic pressure p and acoustic volume flow Û along the instrument's longitudinal coordinate , assuming time-harmonic oscillation with angular frequency ω > 0:

       d p d (iω, ) + Z v (iω, ) Û (iω, ) = 0, d Û d (iω, ) + Y t (iω, ) p(iω, ) = 0. (1)
This formulation of the equations highlights their "transmission line" form [2, (5.132)]. It is also possible to write this model as a second order equation on the pressure only, often attributed to Webster [START_REF] Webster | Acoustical impedance and the theory of horns and of the phonograph[END_REF] although it can be found in much earlier works [START_REF] Bernoulli | Recherches physiques mécaniques et analytiques sur le son et sur le ton des tuyaux d'orgues différemment construits[END_REF][START_REF] Lagrange | Nouvelles recherches sur la nature et la propagation du son (New research on the nature and propagation of sound)[END_REF], as noted in [START_REF] Eisner | Complete Solutions of the "Webster" Horn Equation[END_REF]. When neglecting diffusion phenomena, the lineic impedance Z v and the shunt admittance Y t are purely imaginary:

Z v,lossless (iω, ) = iω ρ 0 S( ) , Y t,lossless (iω, ) = iω S( ) ρ 0 c 2 0 ,
where S( ) is the cross-section area at abscissa , ρ 0 is the static density of air, and c 0 is the sound velocity, which for ideal gases is expressed as γp 0 /ρ 0 with p 0 the static pressure [2, (5.132), (1.98)].

The main limitations of this model are the following: firstly, the model does not take into account non-planar modes of propagation, which limits its validity to tubes with slow longitudinal variation of section, and to wavelengths much longer than the tube diameter; secondly, it does not take into account energy dissipation due to viscous and thermal effects, and therefore predicts inaccurately the frequency, the width and the amplitude of the air column resonances.

Several techniques are available to obtain better predictions in wide tubes or in horns with abrupt section variation [2, §7.6.3.5] including 3D finite elements [START_REF] Rucz | Innovative methods for the sound design of organ pipes[END_REF], multi-modal propagation [START_REF] Pagneux | A study of wave propagation in varying cross-section waveguides by modal decomposition. Part I. Theory and validation[END_REF][START_REF] Guennoc | Improved multimodal formulation of the wave propagation in a 3D waveguide with varying cross-section and curvature[END_REF], and ad hoc methods [START_REF] Fletcher | Acoustic impedance measurements-correction for probe geometry mismatch[END_REF]. The use of a horn equation with a modified effective radius profile has also shown improved results [START_REF] Hélie | Unidimensional models of acoustic propagation in axisymmetric waveguides[END_REF][START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF][START_REF] Eveno | Wave Propagation and Radiation in a Horn: Comparisons Between Models and Measurements[END_REF]. This issue is outside the scope of the present article, which instead focuses on the viscous and thermal phenomena, and in which the waves are assumed to be only planar or spherical.

In order to take into account viscous and thermal effects, which mainly manifest by forming boundary layers near the tube wall, the lineic immittances Z v and Y t are modified.

Z v (iω, ) = 1 1 -K v (iω, ) iωρ 0 S( ) , Y t (iω, ) = 1 + (γ -1) K t (iω, ) iω S( ) ρ 0 c 2 0 . (2)
The coefficient K v describes the perturbation of the lossless lineic impedance due to the viscosity, and K t describes the perturbation of the lossless shunt admittance due to thermal effects. The influence of the boundary layers on the acoustic propagation is contained in the expression of these coefficients.

Note that from the theory of electrical transmission lines, it is possible to interpret the coefficients Z v (iω, ) and Y t (iω, ) as series impedance and shunt admittance per unit length [START_REF] Chaigne | Acoustics of musical instruments[END_REF] 1 .

In that case, the definitions

Γ(iω, ) = Z v (iω, )Y t (iω, ) and Z c (iω, ) = Z v (iω, )/Y t (iω, )
correspond respectively to a local propagation number (purely imaginary in the absence of dissipation) and a local characteristic impedance (real in the absence of dissipation). In general, their real and imaginary parts will be nonzero and depend on frequency and space. System (1) can be solved using the transfer matrix method on intervals where the lineic immittances are constant with respect to space [START_REF] Chaigne | Acoustics of musical instruments[END_REF], or using 1D finite element methods in the general case [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF][START_REF] Chabassier | Open wind instrument design -a python toolbox assisting instrument makers[END_REF].

Other methods exist for computing acoustic propagation with viscous and thermal effects.

The exact solution in a rigid cylinder was derived by Kirchhoff [START_REF] Kirchhoff | Ueber den Einfluss der Wärmeleitung in einem Gase auf die Schallbewegung (On the influence of heat conduction in a gas on the motion of sound)[END_REF]; however, its numerical computation requires an iterative solving of the implicit dispersion relation. The linearized Navier-Stokes equations may directly be solved in arbitrary geometries using 3D finite elements [START_REF] Kampinga | Performance of several viscothermal acoustic finite elements[END_REF]. Disregarding bulk losses leads to a more efficient 3D model called Sequential Linearized Navier-Stokes (SLNS) [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF], with negligible error in the audible range, which is used as a reference 1 Their real part should be nonnegative at all frequencies since they represent a passive system [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF] in the numerical comparisons of section 7. When the boundary layers are thin with respect to the radius of curvature of the walls, equivalent boundary conditions have been found ; this reduces the problem to a Helmholtz equation with effective wall admittance [START_REF] Cremer | On the acoustic boundary layer outside a rigid wall[END_REF][START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF].

Various authors have focused on the special case of infinite right circular cylinders with rigid and isothermal walls, and have derived approximate models in various regimes of frequency and/or pipe size [START_REF] Kirchhoff | Ueber den Einfluss der Wärmeleitung in einem Gase auf die Schallbewegung (On the influence of heat conduction in a gas on the motion of sound)[END_REF][START_REF] Rayleigh | The theory of sound[END_REF][START_REF] Weston | The theory of the propagation of plane sound waves in tubes[END_REF][START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: Transmission line parameter approximations for isothermal and nonisothermal boundary conditions[END_REF][START_REF] Scheichl | On the calculation of the transmission line parameters for long tubes using the method of multiple scales[END_REF], often giving rise to fractional-order derivatives. A review of many such models was done by Tijdeman [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF]. The model initially presented by Zwikker and Kosten [START_REF] Zwikker | Sound Absorbing Materials[END_REF][START_REF] Iberall | Attenuation of oscillatory pressures in instrument lines[END_REF], which shall be recalled in section 6, is valid for most of the frequencies and tube sizes considered in musical acoustics. Tijdeman [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF] calls this model the "low reduced frequency solution", as it is valid as soon as the "reduced frequency" kR = Rω/c 0 is small and the tube is long enough2 (see Table 1 for notations and typical values). Then the transmission line coefficients only depend on the shear wave number3 Sh = R ρ 0 ω/µ. This model, widely used in the musical acoustics community to describe straight circular pipes, is also used for tubes with varying radius [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF].

The wide spread of this debatable practice led us to question whether the Zwikker-Kosten model of viscothermal effects remains reliable for arbitrary variations of radius. Indeed, the influence a nonconstant section would have on the viscothermal effects is unclear. We examine this question in the particular case of cones.

The current work builds upon a previous generalization of the Zwikker-Kosten model to straight tubes with arbitrary cross-sectional shape [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF], which is now frequently used [START_REF] Christensen | Modeling the effects of viscosity and thermal conduction on acoustic propagation in rigid tubes with various cross-sectional shapes[END_REF], especially in the study of porous materials [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF][START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF].

The article is organized as follows. Section 2 recalls the linearization of the Navier-Stokes equations of a compressible fluid for small sinusoidal oscillations. The simplifying assumptions are presented and discussed in section 3. From the simplified equations, the derivation of the corresponding 1D equation is detailed in section 4. The coefficients of this equation are computed in section 5 in the case of cylinders, where the classical Zwikker-Kosten model is recovered; the closed-form solution of system (1) in the whole tube is deduced. Section 6 focuses on the case of circular cones and leads to the main result: a different expression of the loss coefficients. Several discussions ensue, firstly on the integration of the transmission line equations, secondly on the effective loss radius with the proposition of an approximate model. In section 7, impedances calculated using these different models are then compared numerically with a 3D model. Appendix A discusses the passivity of the immittance operators of cylinders and cones in the Laplace domain.

Base equations

The Navier-Stokes equations describing the motion of a perfect gas are linearized for sinusoidal oscillations around a uniform state of air at rest, with no mean flow and small variations of temperature, density and pressure relative to their mean value. In the following, computations are first derived in the frequency domain, classically used in acoustics; the final results are also recast in the complex Laplace domain to analyze causality and passivity in Appendix A.

Let ω > 0 denote a fixed angular frequency. The time-domain fields of pressure P , air velocity v, density ρ and temperature T are assumed to be of the form

P (x, t) = p 0 + Re( P (x)e iωt ), ρ(x, t) = ρ 0 + Re(ρ(x)e iωt ), (3a) v
(x, t) = 0 + Re(v(x)e iωt ), T (x, t) = T 0 + Re( T (x)e iωt ). (3b)
with respect to the 3D spatial variable x and time variable t, where the complex unknowns P (x), v(x), ρ(x), T (x) express the amplitude and relative phase of variations of pressure, velocity, density and temperature around the equilibrium. Static pressure, static density and static temperature are denoted as p 0 , ρ 0 and T 0 , and static flow is assumed to be zero. In all these equations, i denotes the imaginary unit such that i 2 = -1, and Re(•) the real part of a complex number. Under these hypotheses and considering air as a perfect gas, the Harmonic Linearized Navier-Stokes system without source is obtained [START_REF] Bruneau | Fundamentals of acoustics[END_REF][START_REF] Bruneau | General formulation of the dispersion equation in bounded visco-thermal fluid, and application to some simple geometries[END_REF]]

                     iω ρ + ρ 0 ∇ • v = 0, ρ 0 iω v = -∇ P + µ∇ 2 v + ζ + µ 3 ∇(∇ • v), ρ 0 C p iω T = κ∇ 2 T + iω P , P p 0 = ρ ρ 0 + T T 0 , (4a) (4b)
(4c) ( 4d
)
where µ is the shear viscosity, ζ is the bulk viscosity, and κ is the thermal conductivity (see values in table 1). Note that the classical Helmholtz equation of acoustics can be obtained by neglecting viscous and thermal effects (µ = ζ = 0, κ = 0) and eliminating unknowns v and T .

The unknowns are defined on a domain Ω ⊂ R 3 representing a tube; the boundary is composed of a wall Σ, and two end surfaces Γ ± . The wall is assumed to be rigid with a no-slip condition for the gas, and isothermal, namely v = 0 and T = 0 on Σ.

(5)

These "perfect wall" conditions mean that the wall is impermeable and conducts heat much better than air. For porous or flexible materials, or heat insulators, these conditions must be relaxed, for instance using a wall impedance [START_REF] Hélie | Unidimensional models of acoustic propagation in axisymmetric waveguides[END_REF][START_REF] Rienstra | Webster's horn equation revisited[END_REF][START_REF] Boutin | Modèle de propagation acoustique dans un tuyau cylindrique à paroi poreuse (Model of acoustic propagation in a cylindrical pipe with porous wall)[END_REF]. Various boundary conditions are possible on end surfaces Γ ± , for instance a Dirichlet condition on pressure [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF] or coupling with other systems. In the current article, the domain Ω is assumed to be either:

• A straight tube of cross-section Ω c ⊂ R 2 , extending from z -to z + Ω = (x, y, z) ∈ R 3 | (x, y, z) ∈ Ω c × [z -, z + ] . (6)
The boundary of the straight tube is composed of the wall Σ = ∂Ω c × [z -, z + ], and two

planar end surfaces Γ -= Ω c × {z -} and Γ + = Ω c × {z + }.
The longitudinal variable is = z. We define the cross-section at abscissa as Ω = Ω c × { }.

• A cone with apex at the origin, extending from radius r -to r + :

Ω = Ψ(r, θ, φ) ∈ R 3 | (r, θ, φ) ∈ [r -, r + ] × Ω c , (7)
where

Ω c ⊂ [0; π] × [0; 2π)
is the set of angles (θ, φ) contained in the cone, and Ψ is the conversion function from polar to Cartesian coordinates:

Ψ(r, θ, φ) = r sin θ cos φ, sin θ sin φ, cos θ T . ( 8
)
The boundary of the cone is composed of the wall Σ = Ψ [r -, r + ]×∂Ω c , and two spherical

end surfaces Γ -= Ψ {r -} × Ω c and Γ + = Ψ {r + } × Ω c .
The longitudinal variable is = r. The cross-section at abscissa is defined as

Ω := Ψ({ } × Ω c ).
In the special case of circular cones, we will set Ω c = [0, Θ] × [0, 2π); however the methods proposed below can be applied to any smooth shape.

These two types of geometries are represented in figure 1. A parallel is drawn throughout the article between these two cases.

Low reduced frequency approximation

The low reduced frequency approximation consists in simplifying system (4), by separating the propagative behavior which occurs mainly in the longitudinal direction from the dissipative behavior which occurs mainly in the transverse direction. The approximation was initially formulated for circular cylinders [START_REF] Zwikker | Sound Absorbing Materials[END_REF], then generalized to straight tubes with arbitrary cross-section [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF].

This section generalizes it to conical tubes with arbitrary cross-section shape. The simplifications are presented first, then their validity is discussed.

Simplifying assumptions

(A1) The propagating pressure waves are assumed to be purely planar for straight tubes, or purely spherical for cones (no higher-order transverse modes). Formally, this means that the pressure field P can be expressed as a function of one scalar coordinate.

For straight tubes:

P (x) = p(z) ∇ P = d p dz e z . (9)
For cones:

P (x) = p(r), ∇ P = d p dr e r , (10)
where r = ||x|| = x 2 + y 2 + z 2 is the radial component in spherical coordinates, and e r = x/||x|| is the radial unit vector. In the following, symbol is used to denote z or r in each case respectively.

(A2) The Laplacian operator ∇ 2 governing viscous and thermal effects in (4b) & (4c) is approximated by the Laplacian along the cross-section ∇ 2 c :

∇ 2 = ∇ 2 c + ∂ 2 ∂ 2 ≈ ∇ 2 c . (11)
This corresponds to neglecting the longitudinal diffusion process involved in the dissipation phenomena.

(A3) The term involving the gradient-of-divergence of v is set to zero:

ζ + µ 3 ∇(∇ • v) |∇p| . (12)
This corresponds to neglecting losses due to dilation.

After these approximations, system (4) simplifies into:

                     iω ρ + ρ 0 ∇ • v = 0, ρ 0 iωv = - d p d + µ∇ 2 c v, ρ 0 C p iω T = κ∇ 2 c T + iω p, p p 0 = ρ ρ 0 + T T 0 , (13a) (13b) (13c) (13d)
where the pressure now depends only on , and we have isolated the longitudinal velocity

v(x) = v(x) • e . (14)
In section 4, solutions of equation ( 13) are shown to satisfy a one-dimensional partial differential equation.

Discussion of the simplifications (A1)-(A3)

Several characteristic lengths appear in system (4):

• λ = 2πc 0 /ω the wavelength in the absence of losses, where c 0 = γp 0 /ρ 0 is the speed of sound in the absence of losses,

• δ v = µ/ρ 0 ω the characteristic length of viscous effects, i.e. viscous boundary layer thickness,

• δ t = κ/ρ 0 C p ω the thermal boundary layer thickness.

• a characteristic length R of the domain in the transverse direction (e.g. the radius for a cylinder)

The ratios of these lengths are natural dimensionless parameters related to the qualitative behavior of system (4):

• kR = 2πR/λ = Rω/c 0 is called the reduced frequency,

• Sh = R/δ v = R ρ 0 ω/µ is called the shear wave number, and

• Pr = δ 2 v /δ 2 t = µC p /κ is called the Prandtl number.
Together with the heat capacity ratio γ, they fully characterize viscothermal propagation in cylinders. This nondimensionalization allows to compare different models [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF][START_REF] Thibault | Viscothermal models for wind musical instruments[END_REF]. For a given gas, coefficients Pr and γ can be considered constant (see Table 1), whereas kR and Sh depend on both pipe radius and frequency.

Note that for brevity, the dependency of k, λ, Sh, δ v and δ t on the angular frequency ω is omitted except when needed.

The validity of simplifications (A1)-(A3) can also be expressed using kR and Sh. • The assumption of planar or spherical waves (A1) is often adopted when studying wind instruments using a one-dimensional approach [START_REF] Bilbao | Direct simulation of reed wind instruments[END_REF][START_REF] Chaigne | Acoustics of musical instruments[END_REF][START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF]. It has been shown to be asymptotically valid when the tube radius is much smaller than the wavelength [START_REF] Rienstra | Webster's horn equation revisited[END_REF]; or in other words, when the reduced frequency is small:

R λ i.e. kR 1.
As illustrated in figure 2, this assumption is invalidated for wide tubes or at high frequencies.

• One possible justification [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF] of simplification (A2) is based on the comparison of the relative magnitudes of the different terms of (4b) and (4c). The variation of v along direction is expected to be mostly due to acoustic propagation, and thus to vary with wavenumber k; thus the amplitude of term

µ d 2 v d 2 is approximately µ k 2 |v|.
This quantity is negligible compared to the left-hand side of (4b) when µ k 2 ρ 0 ω, i.e. when the viscous boundary layer thickness is much smaller than the wavelength:

δ 2 v k -2 , i.e. (kR) 2 / Sh 2 1. (15)
Similarly, κ d 2 T d 2 is negligible compared to ρ 0 C p iω T when the thermal boundary layer thickness is much smaller than the wavelength.

δ 2 t k -2 i.e. (kR) 2 /( Pr Sh 2 ) 1, (16)
where the Prandtl number Pr is close to unity (see Table 1).

As can be visualized in Figure 2, the ratio δ v /λ and consequently (kR) 2 / Sh 2 are less than 1 × 10 -3 at all audible frequencies, which validates this simplification.

• A justification of assumption (A3) is given in [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF]. The relative magnitudes of the right-hand side terms of (4b) are estimated, by using the conservation of mass, assuming that ∇ρ/ρ 0 is of the same order as ∇p/p 0 (which is true for adiabatic and for isothermal compression), and that ζ is not much bigger than µ:

ζ + µ 3 ∇(∇ • v) = ω ζ + µ 3 |∇ρ/ρ 0 | ≈ ωµ |∇p/p 0 | ≈ (kR) 2 Sh 2 |∇p| .
Using the same estimate as above, at all audible frequencies this term is much smaller than ∇p, and can be neglected. Some authors use this simplification (A3) alone in a 3D finite-element setting, to dissociate the viscous and the thermal effects [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF].

It should be noted that many models of viscothermal acoustics make the assumption of a large shear wave number [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF][START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF][START_REF] Berggren | Acoustic boundary layers as boundary conditions[END_REF], which is invalid for narrow tubes at low frequencies.

Sh 1 i.e. R/δ v 1
This is not the case of the Zwikker-Kosten model, nor of the hereby proposed generalization.

Derivation of the 1D equations

In this section, the simplified system of equations ( 13) is reformulated as a pair of 1D differential equations coupling the pressure and longitudinal velocity. The procedure requires no additional simplification; it is similar to that in [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF] for straight tubes, but extends it to conical tubes with arbitrary cross-sectional shape.

To derive the 1D equations, one must express the velocity and density in terms of the pressure unknown, and relate them through the equation of conservation of mass. For conciseness, the dependency on ω > 0 is omitted in subsections 4.1 and 4.2.

Shape function

The simplified equations for velocity (13b) and temperature (13c) share the same form: a diffusion equation with source term depending on the pressure. A common dimensionless version with unknown ψ(α, x) is given by (recalling the notations from figure 1 and equation ( 11))

     ψ(α, x) + R 2 ( ) α ∇ 2 c ψ(α, x) = 1 inside the domain x ∈ Ω, ψ(α, x) = 0 on the wall x ∈ Σ, (17)
where α ∈ C is a dimensionless parameter balancing propagation and diffusion effects, and where R ≡ R( ) is a characteristic size of the domain in the transverse direction. Parameter α is set to either α v or α t :

α v ( ) = -i Sh 2 = -R 2 ( )ρ 0 iω/µ, α t ( ) = -i Sh 2 Pr = -R 2 ( )ρ 0 iωC p /κ. (18)
Then the velocity and temperature are related to ψ through the following formulas:

v(x) = - 1 ρ 0 iω d p d ( ) ψ α v ( ), x , T (x) = 1 ρ 0 C p p( ) ψ α t ( ), x , (19)
The unknown ψ is closely related to the scalar viscous and thermal fields used by [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF] in the 3D setting. Using the equation of state (13d) to relate ρ to T , and using that p 0 = γ-1 γ C p ρ 0 T 0 , the density can also be expressed in terms of ψ and p:

ρ(x) = ρ 0 p 0 p( ) 1 - γ -1 γ ψ α t ( ), x . (20)
Remark. The choice of introducing R 2 both in the numerator and the denominator (inside α) of equation ( 17) may seem arbitrary. One justification is that this choice makes α dimensionless: in fact we can identify above that α v ( ) = -i Sh 2 ( ) and α t ( ) = -i Pr Sh 2 ( ). Additionally, it suggests to rescale the domain of equation ( 17), using a change of variable x = x/R.

Since equation [START_REF] Guennoc | Improved multimodal formulation of the wave propagation in a 3D waveguide with varying cross-section and curvature[END_REF] involves no derivative in the longitudinal direction, it can be considered on each cross-section Ω independently. To highlight the fact that the overall shape of ψ depends only on α, the problem is brought back to a reference domain Ω, which is planar for straight tubes, or a section of a sphere for cones. In the latter case we impose that R be proportional to the distance to the apex, so that Ω = 1 R Ω . Then parameter α controls the thickness of the boundary layers in the reference domain. The rescaled equation is:

ψ(α, x) = ψ(α, x/R),      ψ(α, x) + 1 α ∇ 2 c ψ(α, x) = 1 on Ω = 1 R Ω , ψ(α, x) = 0 on ∂ Ω. (21)
In particular, the average value of ψ on Ω depends only on the complex number α. We thus introduce:

F (α) = 1 S( ) Ω ψ(α, x) dS, where S( ) = Ω dS. (22)
Function F (α) fully describes how the shape of the cross-section affects viscothermal losses. This shape function F (α) must be calculated for each considered cross-sectional shape, by solving equation [START_REF] Guennoc | Improved multimodal formulation of the wave propagation in a 3D waveguide with varying cross-section and curvature[END_REF]. Analytical expressions are given below for cylinders and for circular cones (see [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF] and ( 52)).

Effective velocity and density

Using equations ( 19) and ( 20), the total flow through one cross-section, defined as Û ( ) = Ω v dS, and the averaged density, defined as ρ ( ) = ( Ω ρ dS)/S( ), can be related to the pressure as

Û ( ) = - S( ) ρ 0 iω 1 -K v ( ) d p d , (23) ρ
( ) = 1 c 2 0 1 + (γ -1)K t ( ) p( ), (24)
where viscous and thermal loss coefficients

K v ( ) = 1-F (-i Sh 2 ( )) and K t ( ) = 1-F (-i Pr Sh 2 ( ))
are introduced. In the limit of lossless propagation (|α| → ∞) these complex coefficients are zero since ψ(α, x) tends to 1 and therefore so does F (α).

Integrating the equation of conservation of mass (13a) on a cross-section Ω yields:

iω Ω ρ dS + ρ 0 Ω ∇ • v dS = 0. ( 25
)
As no flow crosses the boundary Σ, the integral of ∇ • v amounts to the longitudinal variation of the flow:

Ω ∇ • v dS = d Û d . (26)
Thus, equation ( 25) can be reformulated

iωS ρ + ρ 0 d Û d = 0. (27)
This equation relates the mean values of the velocity and of the mass density, and involves only one spatial variable .

Resulting 1D equation

Two relations between the pressure and mean velocity are now available: one is [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF], and the second one is obtained by replacing the density by its expression ( 24) into [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF].

Finally, this leads to the transmission line system (1) for p and Û , where the lineic immittances are of the form given in equation ( 2):

       d p d (iω, ) + 1 1 -K v (iω, ) iωρ 0 S( ) Û (iω, ) = 0, d Û d (iω, ) + 1 + (γ -1) K t (iω, ) iω S( ) ρ 0 c 2 0 p(iω, ) = 0. (28)
with loss coefficients defined by

K v (iω, ) = 1 -F (-i Sh 2 (ω, )) and K t (iω, ) = 1 -F (-i Pr Sh 2 (ω, )) (29)
All the information about the viscothermal losses is contained in coefficients K v and K t ; their computation requires to evaluate the shape function. In the general case, this function can be computed by solving equation ( 17) with numerical methods. Some special cases lead to analytic expressions. Two of them are derived hereafter: in the case of circular cylinders, which is well known, and in the case of circular cones, which had not yet been studied as far as the authors know.

Low reduced frequency solution for right circular cylinders

The procedure of section 4 is now applied to the classical case of right circular cylinders. The calculation of the shape function defined in [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF] leads back to the well-known Zwikker-Kosten model [START_REF] Zwikker | Sound Absorbing Materials[END_REF][START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF][START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF][START_REF] Chaigne | Acoustics of musical instruments[END_REF]. Let R denote the radius of the cylinder. In all this section = z ∈ [z -, z + ] denotes the axial coordinate, and r ∈ [0, R] the radial coordinate in the cylindrical basis.

Shape function of a cylinder

Using the notations of section 3, the cross-section is Ω c = {(x, y) | x 2 + y 2 ≤ R 2 }, and the considered domain is

Ω = (x, y, z) ∈ R 3 | (x, y, z) ∈ Ω c × [z -, z + ] (30)
Due to symmetry, the solution is a function ψ(α, r) only depending on the radial coordinate.

     ψ(α, r) + R 2 α 1 r ∂ ∂r r ∂ ψ ∂r = 1 on Ω, ψ(α, R) = 0. (31)
Using the change of variable ψ(α, r) = ψ(α, rR/ √ α) -1, where r = r √ α/R is a dimensionless coordinate, this equation can be identified with a Bessel equation:

ψ + 1 r ∂ ∂ r r ∂ ψ ∂ r = 0.
Hence, ψ is of the form

ψ(α, r) = 1 + C 1 J 0 √ α r R + C 2 Y 0 √ α r R , (32)
where J 0 and Y 0 are Bessel functions of the first and second kind, and C 1 and C 2 are two complex constants. Since regular solutions are sought, ψ must be continuous at r = 0, so that C 2 = 0.

The other constant is obtained from the boundary condition ψ = 0 at r = R that yields

ψ(α, r) = 1 - J 0 ( √ αr/R) J 0 ( √ α) . ( 33
)
For a given α ∈ C, the average value of ψ over Ω is

ψ = R 0 2πr ψ(r)dr R 0 2πr dr . ( 34
)
Using the Bessel function identity

d dx (xJ 1 (x)) = xJ 0 (x),
allows to derive this average value

F cyl (α) = ψ = 1 - 2J 1 ( √ α) √ αJ 0 ( √ α) . ( 35
)
Function F cyl (α) is the shape function of the cylinder: it expresses the averaged behavior of viscous and thermal diffusion processes on a given cross-section.

Note that due to parity properties of Bessel functions, expression [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF] does not depend on which determination of the square root is chosen. Namely, the power series for J 0 ( √ α) and

2J 1 ( √ α)/ √ α around 0 are respectively [43] J 0 ( √ α) = ∞ m=0 (-1) m (m!) 2 α 4 m , 2J 1 ( √ α)/ √ α = ∞ m=0 (-1) m m!(m + 1)! α 4 m . ( 36
)
Both series have an infinite radius of convergence, thus F cyl (α) is analytic except where α is the square of a zero of J 0 :

α ∈ C \ {j 2 k | (j k ) k∈N zeros of J 0 } Function F (α)
is therefore analytic everywhere except on a discrete set of poles located on the positive part of the real axis. Since α is purely imaginary (either -iωR 2 ρ 0 /µ or -iωR 2 ρC p /κ), function F depends smoothly on the frequency ω > 0.

Transmission line equations for the cylinder

The transmission line equation is rewritten using the shape function of the cylinder [START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF]:

         1 1 -K v,cyl (iω) iωρ 0 πR 2 Û (z) + d p dz (z) = 0, 1 + (γ -1)K t,cyl (iω) iω πR 2 ρ 0 c 2 0 p(z) + d Û dz (z) = 0, (37)
where the loss coefficients do not depend on z and are given by

∀ω > 0, K v,cyl (iω) = J(R -iωρ 0 /µ) = J(i 3/2 Sh(ω)), (38)
K t,cyl (iω) = J(R -iωρ 0 C p /κ) = J(i 3/2 √ Pr Sh(ω)), (39)
with

J(q) = 2J 1 (q) q J 0 (q) . (40)
System [START_REF] Christensen | Modeling the effects of viscosity and thermal conduction on acoustic propagation in rigid tubes with various cross-sectional shapes[END_REF] coincides with the model initially derived by Zwikker and Kosten [START_REF] Zwikker | Sound Absorbing Materials[END_REF][START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF][START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary cross-sectional shape[END_REF]. In its "transmission line" form, it has been used by various authors [2, (5.132)], [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF][START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Ernoult | Transfer matrix of a truncated cone with viscothermal losses: application of the WKB method[END_REF][START_REF] Boutin | Modèle de propagation acoustique dans un tuyau cylindrique à paroi poreuse (Model of acoustic propagation in a cylindrical pipe with porous wall)[END_REF]. The extension to the Laplace domain of the lossy immittance operators appearing in equation ( 37) is discussed in Appendix A.1.

Propagation constant and characteristic impedance

Since coefficients do not depend on the axial coordinate z = , the solutions of (37) are given by the classical solutions of the transmission line equations:

p(z) = p+ e -Γz + p-e Γz , Û (z) = 1/Z c p+ e -Γz -p-e Γz , (41)
where the complex propagation constant and characteristic impedance are given by

Γ = iω c 0 1 + (γ -1)K t,cyl (iω) 1 -K v,cyl (iω) , Z c = ρ 0 c 0 πR 2 1 -K v,cyl (iω) 1 + (γ -1)K t,cyl (iω) -1/2 (42)
and p+ , p-are the complex amplitudes of the "forward" and "backward" wave respectively. Note that the loss coefficients K v and K t depend on frequency, resulting in wave dispersion as Γ depends non-linearly on ω. Moreover, Γ has a positive real part accounting for the wave attenuation with respect to distance.

Asymptotic expansions of K v and K t were proposed by Keefe [START_REF] Keefe | Acoustical wave propagation in cylindrical ducts: Transmission line parameter approximations for isothermal and nonisothermal boundary conditions[END_REF] for large shear wave number Sh, i.e. when the boundary layers are thin compared to the tube's radius. Various expressions of the dispersion relation due to viscothermal effects in a cylinder have been proposed in the past. A review by Tijdeman [START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF] showed that most of them are approximations of the Zwikker-Kosten model.

Low reduced frequency solution for right circular cones

We now follow exactly the same steps as in section 5, adapted to the case of a cone with circular base. The solution ψ of equation ( 17) is calculated for this geometry; its average yields the shape function F (α) defined in [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF], resulting in a new expression of loss coefficients. Let Θ denote the half opening angle of the cone. In all this section = r ∈ [r -, r + ] denotes the distance of one considered point from the apex, and θ ∈ [0, Θ] the angle measured from the cone axis, in spherical coordinates (see Figure 3). Using the notations of section 2, the set of angles contained in the cone is Note that equation ( 17) makes use of a characteristic transverse length. To reflect the widen-ing of the cone with the distance to the apex, we set

Ω c = [0, Θ] × [0, 2π).
R(x) = c Θ r, (43)
where the constant c Θ will be chosen later. It should however be close to Θ for small cone angles

(c Θ = Θ + O(Θ 2 )).
Possible choices for c Θ are discussed in section 6.4.

Shape function of a right circular cone

The shape function F (α) of a cone must be determined by computing the solution ψ(α, x) of equation [START_REF] Guennoc | Improved multimodal formulation of the wave propagation in a 3D waveguide with varying cross-section and curvature[END_REF]. Due to symmetry of the problem around the main axis of the cone, the equation reduces to:

     ψ(α, r, θ) + c 2 Θ α sin θ ∂ ∂θ sin θ ∂ψ(α, r, θ) ∂θ = 1 for 0 ≤ θ ≤ Θ, ψ(α, r, Θ) = 0. ( 44
)
This is an ordinary differential equation in variable θ. The dependence of ψ on α and r is made implicit in the following. Equation ( 44) can be identified as a noninteger-order equation of spherical harmonics:

η(η + 1)ψ + 1 sin θ ∂ ∂θ sin θ ∂ψ ∂θ = η(η + 1), (45)
where η(α, r) ∈ C is such that

η(η + 1) = α c 2 Θ , i.e. η = ± α c 2 Θ + 1 4 - 1 2 . ( 46
)
A classical way of solving this equation is to perform a change of variable; with q := cos θ, the equation ( 45) is rewritten as

d dq (1 -q 2 ) dψ dq + η(η + 1)ψ = η(η + 1). ( 47
)
which is the Legendre differential equation, a special case of hypergeometric equation, for which the solutions are widely known [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. The solutions ψ(θ) have the form:

ψ(θ) = 1 + C 1 P η (cos θ) + C 2 Q η (cos θ), (48)
where P η and Q η denote the Legendre functions of noninteger degree of the first and second kind.

Note that due to properties of the Legendre functions, this expression does not depend on which root is chosen in [START_REF] Thibault | Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments[END_REF]. As in section 5.1, since regular solutions are sought, ψ must be continuous at θ = 0, so that C 2 = 0. The other constant is obtained from the boundary condition ψ(Θ) = 0 which yields

ψ(θ) = 1 - P η (cos θ) P η (cos Θ) . ( 49
)
The shape function can be computed as

F Θ (α) = ψ = Θ 0 r 2 sin(θ) ψ(θ) dθ Θ 0 r 2 sin(θ) dθ . ( 50
)
Legendre function P η satisfies the following useful identity:

Θ 0 sin(θ)P η (cos θ) dθ = sin 2 Θ η(η + 1) P η (cos Θ), (51)
where P η (z) is the derivative of P η (z) with respect to its argument z. Therefore

F Θ (α) = 1 - sin 2 Θ 1 -cos Θ P η (cos Θ) η(η + 1)P η (cos Θ) , where η η + 1 = α c 2 Θ . (52)
This shape function F Θ (α) is parametrized by the cone angle Θ. A calculation can show that in the limit of a small cone angle, the cylindrical shape function is recovered4 .

Transmission line equations of the cone

Pressure and velocity are thus solutions to equations ( 1) and ( 2), where the loss coefficients are calculated using shape function F Θ [START_REF] Chabassier | Direct computation of modal parameters for musical wind instruments[END_REF]. Dependency on position r and frequency ω is highlighted for convenience:

         1 1 -K v,Θ (r, iω) iωρ 0 S(r) Û (r) + d p d (r) = 0, 1 + (γ -1)K t,Θ (r, iω) iωS(r) ρ 0 c 2 0 p(r) + d Û d (r) = 0. for r ∈ [r -, r + ] (53)
The loss coefficients of the cone are given by:

K v,Θ (r, iω) = P η v (r, iω), cos Θ , η v (r, iω) = - 1 2 ± 1 4 -iωr 2 ρ 0 µ (54) K t,Θ (r, iω) = P η t (r, iω), cos Θ , η t (r, iω) = - 1 2 ± 1 4 -iωr 2 ρ 0 C p κ , (55)
where

P(ν, q) = (1 + q) P ν (q) ν(ν + 1)P ν (q) . (56)
Once again, the values of the Legendre function P η and its derivative do not depend on which determination of the square root is chosen for η [43, §8.2.1]. Although the expression of the loss coefficients is more complicated than previously, the structure of the equation relating p and Û is exactly the same as in the cylindrical case. Note that the loss coefficients depend both on frequency and space, similarly to when authors apply the Zwikker-Kosten model to an instrument with a varying radius [START_REF] Bilbao | Direct simulation of reed wind instruments[END_REF][START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF]. The main difference between the classical formulas and these new coefficients is that the Bessel functions of equation ( 40) have become Legendre functions in equation ( 56). The extension of these new operators to the Laplace domain, and the passivity of the immittances calculated using them in equation ( 2), are discussed in Appendix A.2.

Solving the transmission line equations

Whereas in the case of straight tube a closed form solution to the transmission line equations ( 37) is given as the sum of "forward" and "backward" waves by equation ( 41), this breaks down in the case of conical tubes. Indeed, due to the spatial dependence of the coefficients, there is no exact closed-form solution known for the transmission line equations (53). To circumvent this difficulty, several authors use constant "averaged" loss coefficients, both in the transfer matrix method [2, §7.

K v (r, iω) ≈ Kv (iω) = 1 r + -r - r+ r- K v (r, iω) dr (57)
or it can be set to Kv (iω) = K v (r * , iω) where r * is some radius between r -and r + , for instance the arithmetic mean of the radii of the two ends [START_REF] Caussé | Input impedance of brass musical instruments -comparison between experiment and numerical models[END_REF]. After this approximation, a closed-form solution is available [START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF]. The approximation error can be reduced by subdividing each cone into shorter cones, thus making a piecewise constant approximation of K v , but this method converges slowly [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF]. Other approximations have been proposed for solving the propagation equation with nonconstant loss coefficients, using a perturbation method [START_REF] Nederveen | Acoustical aspects of woodwind instruments[END_REF] or the WKB method [START_REF] Ernoult | Transfer matrix of a truncated cone with viscothermal losses: application of the WKB method[END_REF].

In the current work, no further approximation is necessary. A high-order 1D finite element method is used for solving system (53) with arbitrary accuracy [START_REF] Tournemenne | A comparison of a one-dimensional finite element method and the transfer matrix method for the computation of wind music instrument impedance[END_REF][START_REF] Chabassier | Open wind instrument design -a python toolbox assisting instrument makers[END_REF].

A second difficulty is related to time-domain integration. As in the Zwikker-Kosten model, the dependency of K v and K t on frequency is highly nonlinear, leading to complicated convolution kernels in the time domain. Methods have been devised to deal with this issue using auxiliary variables [START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF][START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Thibault | Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments[END_REF], which may be accommodated to use the newly proposed coefficients of system (53).

On the effective loss radius

Several authors choose to use the formulas for the cylinder coefficients, given by equations [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF], even in the case of conical tubes. This requires the choice of an effective "tube radius" R to be used in the formula for K v,cyl , depending on the cone angle Θ and the distance to the apex r. As we are considering spherical waves rather than planar waves, several definitions of this "effective radius" are possible and lead to slightly different results.

(ALT) Most authors implicitly choose the radius of the circle bounding the considered spherical cap [START_REF] Caussé | Input impedance of brass musical instruments -comparison between experiment and numerical models[END_REF][START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Schmutzhard | Parameter optimisation of a viscothermal time-domain model for wind instruments[END_REF]. This corresponds to setting c Θ,ALT = sin Θ and R ALT = r sin Θ. The length R ALT is called the altitude hereafter, as it is the length of the segment which joins a boundary point to its orthogonal projection on the main axis.

(ARC) Another natural choice is c Θ,ARC = Θ: this would correspond to choosing R ARC = Θr, which can be interpreted geometrically as the length of the arc traced from the axis to a point on the edge of the spherical cap.

(HR) A third possible choice is to ensure that the surface-to-perimeter ratio of the cylinder's cross section is the same as that of the spherical cap. We denote R HR the hydraulic radius of the considered spherical cap, given by the following expressions:

R HR (x) = 2 area perimeter = c Θ,HR r, where c Θ,HR = 2(1 -cos Θ) sin Θ . ( 58
)
This choice is inspired from the modeling of porous media [START_REF] Allard | Propagation of sound in porous media: modelling sound absorbing materials 2e[END_REF], and motivated by the fact that boundary layers are thinning as frequency grows, making losses mainly depend on the ratio between wall surface and air volume. Note that the asymptote c Θ,HR ∼ Θ→0 Θ remains valid.

For each convention conv in {ALT, ARC, HR}, we define the corresponding loss coefficient

K v,conv : ∀ω > 0, K v,conv (iω) = J(R conv -iωρ 0 /µ), (59)
and similarly K t,conv (iω) = J(R conv -iωρ 0 C p /κ), where the function J is defined in equation [START_REF] Bruneau | General formulation of the dispersion equation in bounded visco-thermal fluid, and application to some simple geometries[END_REF]. Each of these K v,conv (iω) could serve as an approximation of K v,Θ (iω); in order to determine which of these conventions leads to the least error, a numerical comparison is conducted.

For the sake of the example, the values r = 2 mm and Θ ∈ {π/3, π/10} are chosen. These exaggerated opening angles are used to maximize the discrepancies between the models. Note that although a special case is presented for clarity, the conclusions shall remain valid for other values of r, as well as for the thermal loss coefficient K t , which has a similar expression 5 .

The viscous loss coefficients for cylinders [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] and their counterpart for cones (54) are compared numerically as functions of frequency, for the conventions (ALT) and (HR). The following discussion is also valid for the behavior of thermal loss coefficient K t , which is given by the same formula up to a parameter scaling (replacing Sh with √ Pr Sh).

Legendre function P λ and its derivative for complex order λ are computed using Python library mpmath [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18[END_REF], through the use of the Gaussian hypergeometric function 2 F 1 :

P λ (z) = 2 F 1 -λ, λ + 1; 1; 1 -z 2 , for |1 -z| < 2, (60)
P λ (z) = λ(λ + 1) 2 2 F 1 -λ + 1, λ + 2; 2; 1 -z 2 , for |1 -z| < 2. ( 61
)
where the Gaussian hypergeometric function is defined as :

2 F 1 (a, b; c; z) = ∞ n=0 a(a + 1) . . . (a + n -1) b(b + 1) . . . (b + n -1) c(c + 1) . . . (c + n -1) z n n! . (62)
The viscous loss coefficient K v,Θ (iω), given by ( 54) is compared numerically with the different possible approximations given by equation (59) for conventions (ALT), (ARC) and (HR) in Figure 4.

In all conventions, at low frequencies the coefficients tend to the real number 1. At high frequencies, all the approximations K v,conv (iω) tend to zero with an asymptotic equivalent in

1-i Rconv 2µ
ρ0ω . Note that these asymptotes differ by a constant factor depending on the chosen convention. Asymptotically, coefficient K v,Θ (iω) follows the (HR) convention at high frequencies, whereas the relative error of the two other conventions remains high. One should therefore avoid approximating K v,Θ (iω) with K v,ALT (iω) as it leads to a slight overestimation of its magnitude, especially at high frequencies. The common use in cones of the Bessel-function-based coefficients should be expected to slightly overestimate the amount of losses. This error can reach 30% for the extreme cone angle Θ = π/3 considered here, although for realistic values Θ ≤ π/10, the error appears to be less than 3%. On the other hand, the relative error between K v,Θ and K v,HR is less than 0.4% in all the cases concerned by musical acoustics (Θ ≤ π/10, r ≥ 2 mm and f ≥ 20 Hz on Figure 4). Moreover, this error decreases at high frequencies whereas the others do not. Therefore the hydraulic radius can be interpreted as a relevant effective radius for viscothermal effects.

Existing algorithms which calculate viscothermal losses using the classical formula [START_REF] Champoux | Dynamic tortuosity and bulk modulus in air-saturated porous media[END_REF] can be modified to take into account conicity by using the hydraulic radius instead of the "altitude length" radius.

R HR = c Θ,HR c Θ,ALT R ALT = 2 1 + cos Θ R ALT , K v,Θ (r, iω) ≈ K v,HR (iω) = J(R HR -iωρ 0 /µ) (63)
This approximation alleviates the need to compute hypergeometric functions, which are more rarely available and more computationally demanding than Bessel functions, while recovering the correct high-frequency asymptotes. of the cone, with a maximal error value of 3.365%, whereas the maximal relative error is below 0.055% for both models SH and ZK-HR. No significant divergence of the 3D and 1D results is observed after the cutoff frequency; this can be attributed to the chosen computational setting which imposes wavefronts to be almost spherical, forcing the validity of assumption (A1).

It is possible to interpret the sign of the error in Figure 6b. For the Zwikker-Kosten model, the error on the amplitude is positive for frequencies at which the amplitude is increasing, and negative where it is decreasing: thus one can expect that ZK slightly underestimates the resonance and antiresonance frequencies. On the other hand, the error on the phase is negative where the phase is positive, and conversely; thus one can expect that ZK slightly underestimates the quality factor of the resonances. This is consistent with the conclusions of section 6.4, which indicated that ZK slightly overestimates the amount of losses.

A calculation of the resonance frequencies and quality factors of the 1D models was performed as described in [START_REF] Chabassier | Direct computation of modal parameters for musical wind instruments[END_REF]. The frequency of the first resonance peak was found to be 2165.07 Hz for model ZK, or 2165.18 Hz for model ZK-HR; a relative difference of 0.005%. The quality factor was more significantly modified: 217.5 for model ZK, and 222.8 for model ZK-HR; a relative difference of 2.3%. Therefore the modification of the loss coefficients proposed in this article appears to increase the quality factor, but to leave the frequency almost unaffected. The fact that the quality factor is lower for ZK is consistent with the interpretation of figure 6b above.

In all these calculations, it is key to use the same definition of the impedance for all models, i.e. equation (64) in which the entry surface Γ -is a spherical cap and not a plane, as the area difference between the two leads to a difference of several percent on the impedance.

Similar results were obtained on a wider frequency range and for a different cone, thus we expect these results to be valid in a more general setting.

Note however that model ZK is already an excellent approximation of the 3D phenomena, and that in many cases there is no need to take into account the effect of the conicity on the loss coefficients. Indeed, other aspects of the acoustic propagation not covered in this paper, such as the curvature of the wavefronts or the radiation impedance, can induce much higher discrepancies on the impedance [START_REF] Eveno | Wave Propagation and Radiation in a Horn: Comparisons Between Models and Measurements[END_REF]. As an example, computations on the complete Vox Humana pipe have been performed with the assumptions of "spherical waves" = r and "planar waves" = z (this amounts to a geometric length correction L r = L z / cos Θ [START_REF] Hélie | Unidimensional models of acoustic propagation in axisymmetric waveguides[END_REF]), and a relative difference of 10% was obtained between the first resonance frequencies (local maximum of |Z in |) of each model.

Also note that the computation of model SH requires the evaluation of Legendre functions of complex order, which are costly to evaluate numerically using series summation. Indeed, we were unable to obtain results for model SH for frequencies 3950 Hz and above due to very slow convergence of the series. This further justifies to use model ZK-HR instead.

The conclusion of this numerical comparison is that, in cones, the main difference between the 3D Sequential Linearized Navier-Stokes model and the 1D Zwikker-Kosten model indeed amounts to the modification of the loss coefficients pointed out in this work. Moreover, the approximation proposed in equation ( 63) leads to an improvement of the results by more than an order of magnitude, with no additional cost compared to the direct use of the Zwikker-Kosten model.

Conclusions

This article has shown that the low reduced frequency assumption, used to reduce viscothermal acoustic propagation in a straight tube to a 1D transmission line model relating pressure and flow, can also be employed in conical tubes, leading to similar results. Moreover, it has been established that in both cases, the transmission line coefficients can be expressed using the solution of a heat equation on the cross-section surface -whether this surface is planar, or a portion of a sphere. This makes it theoretically possible to study acoustic losses in tubes or cones, with any cross-sectional shape.

Closed-form solutions of this heat equation have been found for the simplest geometries (right circular cylinders and right circular cones). Numerical comparisons show that the corresponding where P(ν, q) is defined by equation (56). A complementary theoretical study would be necessary to determine the domain on which this analytic extension remains valid; we conjecture that, like in Appendix A.1, the only singularities are a discrete, countably infinite set of strictly negative real poles. This property appears to hold numerically. Like in Appendix A.1, we thus expect (without proof) that K v,Θ (s) and K t,Θ (s) are the Laplace transforms of causal decreasing signals, corresponding to causal and stable operators [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF]VI,Th. 5].

This further allows to extend the lineic impedance Z v and shunt admittance Y t , defined by

(2), to the Laplace domain. It can then be checked that these operators are plausible from a physical standpoint; more precisely, the lineic immittances should be causal and passive operators [START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF]. Figure A.7 plots the analytic extension of these lineic immittances to the complex plane, as a function of the Laplace variable s. Then we observe that, in the left-hand side of the complex plane they display poles (the bright spots in the magnitude plot) and zeros (the dark spots).

On the phase plot, the curve on which the real part of the lineic immittances is equal to zero is highlighted, showing that in the right-hand side of the complex plane, their real part remains Thanks to these properties, the transmission line equations give rise to a model which dissipates energy, and would thus be suitable for use in time-domain simulations using an appropriate finite-dimensional representation [START_REF] Hélie | Diffusive representations for the analysis and simulation of flared acoustic pipes with visco-thermal losses[END_REF][START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF][START_REF] Monteghetti | Design of broadband time-domain impedance boundary conditions using the oscillatory-diffusive representation of acoustical models[END_REF][START_REF] Thibault | Dissipative time-domain one-dimensional model for viscothermal acoustic propagation in wind instruments[END_REF].

  × 10 -5 kg m -1 s -1 Bulk viscosity ζ 1.3 × 10 -5 kg m -1 s -1

2 rad Frequency f 20 to 20 × 10 3 Angular frequency ω 2πf 10 2 to 10 5 rad s - 1 Radius of the pipe R 10 -3 to 10 -1 m Wavelength λ c0/f 2 × 10 -2 to 2 × 10 1 m 3 m

 203110210213 Hz Viscous boundary layer thickness δv µ/(ρ0ω) 10 -5 to 10 -Thermal boundary layer thickness δt κ/(ρ0Cpω) 10 -5 to 10 -3 m Reduced frequency kR Rω/c0 3 × 10 -4 to 30 (no unit) Shear wave number Sh R ρ0ω/µ 3 to 8000 (no unit) Prandtl number Pr µCp/κ 0.71 (no unit) Heat capacity ratio γ Cp/Cv 1.402 (no unit)

Figure 1 :

 1 Figure 1: Sketches of the two considered types of 3D domains, with the notations used in each. Left: a straight tube with arbitrary cross-section ranging from z -to z + . Right: a cone with arbitrary cross-section ranging from radius r -to r + .

Figure 2 :

 2 Figure 2: Range of validity of the low reduced frequency model (grayed), compared to the range of tube radii and frequencies of interest in musical acoustics (rectangle). Iso-lines of the ratios of radius R, wavelength λ and boundary layer thickness δv are shown. The low reduced frequency model assumes R λ and δv λ, but does not assume R δv.

Figure 3 :

 3 Figure 3: Sketch of the 3D domain considered in section 6, with the spherical coordinates used.

4 . 5 ]

 45 , and in time-domain methods[19, §5.2]. The constant value can be chosen to be the mean value of the coefficients[START_REF] Hélie | One-dimensional acoustic models of horns and comparison with measurements[END_REF]

Figure 4 :

 4 Figure 4: Comparison of the loss coefficient of a cone of half-angle Θ = π/3 (top) or π/10 (bottom) at distance r = 2 mm from the apex, with the loss coefficients of cylinders with different radii. The radii of the cylinders are chosen to coincide with the transverse characteristic length of the cone using different conventions (ALT), (ARC) and (HR). The magnitude (left) and phase (center) of each of these models, and the magnitude of the relative error of each approximation (right) are shown as functions of the frequency.

Figure 5 :

 5 Figure 5: Geometry of the considered cone (right). The solid lines represent the instrument wall, and the dotted lines represent the input and exit surfaces.

Figure A. 7 :

 7 Figure A.7: Amplitude (left) and phase (right) of the analytic extension to the complex plane of the lineic immittances Zv (top) and Yt (bottom) defined by equation (2). The loss coefficients (54) have been calculated for a cone with angle Θ = π/4. The calculation uses unitary coefficients R = Pr = ρ 0 = µ = 1. Black curves on the phase plots indicate the locus of points where the real part of the lineic immittances is zero.

  positive: ∀s s.t. Re(s) ≥ 0, Re(Z v (s)) ≥ 0 and Re(Y t (s)) ≥ 0. (A.5)Moreover, both operators are symmetric by complex conjugation.∀s ∈ C, Z v (s) = Z v (s) and Y t (s) = Y t (s) (A.6)

  Two dimensionless complex numbers K v (iω, ) and K t (iω, ) are classically introduced [21, eqns. (3.115-3.117)], termed loss coefficients throughout the paper, related to the viscous and thermal effects respectively.

Although the reduced frequency is denoted k in Tijdeman's article[START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF], the choice has been made here to conform to the literature of acoustics in which k designates the wave number.

This notation and naming comes from[START_REF] Tijdeman | On the propagation of sound waves in cylindrical tubes[END_REF] ; some authors prefer to introduce the two wavenumbers of the viscous and thermal diffusion waves[START_REF] Chaigne | Acoustics of musical instruments[END_REF] [21, (2.85,2.87)], which have the unit of an inverse length.

More precisely, for all α ∈ C such that J 0 (α) = 0, one has lim Θ→0 F Θ (α) = F cyl (α). The proof is not included, but can be summarized as follows: write F cyl (α) and F Θ (α) using ratios of functions; express the numerators and denominators as series of functions; prove pointwise convergence of each term of each series as Θ → 0; for a given α, prove normal convergence of the series for Θ small enough; conclude that the numerator and denominator converge to their expected limit, and thus that F Θ (α) converges to F cyl (α), for all α such that the denominator is nonzero.

The loss coefficients K v,cyl , K t,cyl , K v,Θ and K v,Θ can in fact all be expressed in terms of the dimensionless quantities Sh, Pr and c Θ , where the shear wave number is proportional to the square root of the frequency, the Prandtl number can be assumed to be constant, and c Θ depends only on the cone angle and the chosen convention.

Numerical comparison

In order to assess the magnitude of the correction considered in this article, impedances of cones are computed using different models. We solve the transmission line equations (1) with = r as the abscissa 6 , with lineic immittances given by [START_REF] Chaigne | Acoustics of musical instruments[END_REF], and several possible expressions of the loss coefficients: ZK. (Zwikker-Kosten) the loss coefficients are calculated from equation (59) using convention (ALT); this corresponds to the viscothermal model many authors use [START_REF] Chaigne | Acoustics of musical instruments[END_REF][START_REF] Bilbao | Passive time-domain numerical models of viscothermal wave propagation in acoustic tubes of variable cross section[END_REF].

SH. (Spherical Harmonics) the loss coefficients are those derived in section 6, given by equations (54-55); ZK-HR. (Zwikker-Kosten with Hydraulic Radius) the loss coefficients are calculated using the equation (63) based on the hydraulic radius proposed in section 6.4.

As a reference, we use the 3D Sequential Linearized Navier-Stokes (SLNS) model of viscothermal acoustics described in [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF].

In [START_REF] Rucz | Modeling resonators of reed organ pipes[END_REF], a Vox Humana organ pipe is described in three parts, for which the dimensions are given in Table 2. It is known that the assumption of single-mode propagation (A1) is the largest source of error in the Zwikker-Kosten model, and numerous authors have studied ways to remedy it [2, 7.6.3.5]. As this issue is out of scope of the current article, the chosen computational setting should be compatible with spherical wavefronts, in order to only observe the differences which are due to the transmission line coefficients themselves. This is why the predictions of the different models are compared on a single cone (mid-part of the Vox Humana pipe), limited by two spherical caps. The corresponding geometry is represented in Figure 5.

Using the same notations as in section 2, a constant-pressure condition P (x) = P0 is applied on the first end surface Γ -, and a simplified "free air" condition P = 0 on the second Γ + . The input impedance of the cone is then computed as the ratio between the applied pressure and the resulting total flow through the input surface. 6 Note that choosing the abscissa = z in a cone is also possible, which leads the 1D domain to be shorter length and significantly modifies the result [START_REF] Eveno | Wave Propagation and Radiation in a Horn: Comparisons Between Models and Measurements[END_REF]. Since the goal of this comparison is to perturb only the viscothermal coefficients, we prefer to use abscissa = r in all this section. The values of the physical coefficients are calculated for air at 20 °C, as given in Table 1.

Section

For the 1D models, the calculations are performed using high-order 1D finite elements with the Python library Openwind [START_REF] Chabassier | Open wind instrument design -a python toolbox assisting instrument makers[END_REF], by applying a unitary flow at the tube entrance, and calculating Z in as in equation ( 64). The Python code for the 1D computations is made available online along with the library. The reference 3D solution is calculated with an axisymmetric formulation, using high-order finite elements on a 2D mesh with the software Montjoie [START_REF] Duruflé | Montjoie software[END_REF]; further implementation details are given in [START_REF] Thibault | Viscothermal models for wind musical instruments[END_REF]. The frequency ranges from 2000 Hz to 4500 Hz. This frequency range is chosen as it encompasses the first resonance frequency (2150 Hz) and extends beyond the first anti-resonance frequency of the cone (2625 Hz). A frequency step of 25 Hz is chosen to limit the computational effort. The cutoff frequency of the first non-planar mode of a cylinder with the same radius R 1 = 24.60 mm is f c = 4088 Hz (or 8500 Hz for the first axisymmetric non-planar mode); beyond that frequency, the presence of higher-order modes becomes very likely in practice [2, eq. (7.147)], hence assumption (A1) is invalid and none of the 1D models is usable. In total, the 3D calculations run in about 1 hour on a laptop CPU (Intel Core i7-3687U, 2.10 GHz), whereas the calculations for model SH take about 1 minute, most of which is spent calculating the hypergeometric series (62), and the calculations for models ZK and ZK-HR take around 0.1 seconds. Due to the small relative differences between the results, all plots would overlap if represented on the same figure. Therefore the magnitude and the phase of the impedance obtained with the 3D model are represented in Figure 6a, and the signed relative error of each 1D model is represented in Figure 6b using a symmetrical logarithmic scale.

It can be observed that all 1D models give an accurate prediction of the impedance, with a mean relative error less than 0.204% over the considered frequency range. Additionally, the impedances computed with the two proposed models SK and ZK-HR are extremely close to one another, with a relative difference below 0.005% at all frequencies, confirming the validity of approximation (63). The mean error is 0.204% for model ZK, 0.023% for SH, and 0.031% for ZK-HR. Model ZK appears to be making significant error near the resonances and the antiresonance transmission line coefficients are different depending on the taper of the tube 7 . It is therefore interesting to find how to take into account this variation in the simplest possible way; we propose to use the same formulas as in cylinders, but to replace the altitude radius by the hydraulic radius in the argument of the Bessel functions. This manipulation should be easy to apply to any program which computes lossy acoustic propagation, and restores the correct high-frequency behavior of the loss coefficients; the hydraulic radius can therefore be interpreted as an effective loss radius.

In our numerical experiments, the corrected models significantly outperform the classical model when used to compute impedance curves. Although the "classical" Zwikker-Kosten model is already extremely accurate (which validates the common assumption that the viscothermal effects in a cone can be approximated by those in a cylinder), the correction of the loss coefficients proposed in this paper leads to a reduction of the error on the impedance by more than an order of magnitude.

One limitation of this numerical comparison is the use of the Sequential Linearized Navier-Stokes model as the 3D reference, since it also performs some simplifying assumptions [START_REF] Kampinga | An efficient finite element model for viscothermal acoustics[END_REF] which may affect the computed impedance. To better ascertain the quality of our models, it would be desirable to perform similar comparisons using a finite element implementation of the full Linearized Navier-Stokes system.

A possible extension of this work could be the study of viscothermal acoustics through the lens of the multi-modal method [START_REF] Pagneux | A study of wave propagation in varying cross-section waveguides by modal decomposition. Part I. Theory and validation[END_REF][START_REF] Guennoc | Improved multimodal formulation of the wave propagation in a 3D waveguide with varying cross-section and curvature[END_REF], which computes the acoustic field by splitting it into a sum of transverse modes. Although it has been developed for lossless acoustic propagation, this method may be able to reveal the higher-order influence of the boundary layers on the acoustic propagation. Moreover, such a method would be much more adaptable to realistic musical instruments with flare, possible side holes, and radiating into the external air.

Appendix A. Loss coefficients in the Laplace domain

Since the system is lossy, it would be sensible to consider whether the loss coefficients can be extended to the Laplace domain, to represent linear time-invariant (LTI) operators. The discussion is done separately for the coefficients K v,cyl (iω) and K t,cyl (iω) obtained for a circular cylinder and given by equations ( 38)- [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF], and for the new formulas of K v,Θ (r, iω) and K t,Θ (r, iω)

7 with the usual convention to use the altitude radius in the Bessel function: this is mainly due to the difference in perimeter-to-surface ratio between a circle and a spherical cap with the same external radius obtained in the case of the circular cone and given by equations ( 54)-(55).

Appendix A.1. Cylinder coefficients in the Laplace domain

The calculations leading from system (13) to the transmission line equations ( 37) can be generalized by replacing iω by the complex Laplace variable s, and remain valid except for particular values of s. In particular, the solution of the transverse problem given by equation [START_REF] Scheichl | On the calculation of the transmission line parameters for long tubes using the method of multiple scales[END_REF] becomes invalid only when the denominator J 0 ( √ α) is zero. Away from these poles, function F cyl (α) remains well defined by the same formula. The loss coefficients are therefore also defined for all s ∈ C except on a discrete set of poles, which are all located on the negative real axis, by:

These expressions are symmetric with respect to complex conjugation, and analytic on the half-

The fact that this domain contains the right half-plane {Re(s) ≥ 0}, can be related to the causality of the corresponding operators. Moreover, since it contains a half-plane of the form {Re(s) ≥ -ε} (ε > 0), this domain is the region of convergence of the Laplace transform of a stable operator [START_REF] Schwartz | Mathematics for the Physical Sciences[END_REF]. One can thus expect that K v,cyl and K t,cyl correspond to the Laplace transform of causal, decreasing, real signals 8 .

Appendix A.2. Cone coefficients in the Laplace domain

As in Appendix A.1, note that all the calculations of section 6.1 can be generalized to the Laplace domain, by replacing iω with the complex variable s. Indeed, the solution (49) of equation ( 44) is valid not only for pure imaginary α, but for any value of α such that the denominator P η(α) (cos Θ) is nonzero. The coefficients K v,Θ (r, iω), which have been computed for sinusoidal oscillations only, can therefore be extended analytically to operators in the Laplace domain using the same formulas.