Analysis of the Effects of the Viscous Thermal Losses in the Flute Musical Instruments

: This article presents the third part of a larger project whose final objective is to study and analyse the effects of viscous thermal losses in a flute wind musical instrument. After imple-menting the test bench in the first phase and modelling and validating the dynamic behaviour of the simulator, based on the previously implemented test bench (without considering the losses in the system) in the second phase, this third phase deals with the study of the viscous thermal losses that will be generated within the resonator of the flute. These losses are mainly due to the friction of the air inside the resonator with its boundaries and the changes of the temperature within this medium. They are mainly affected by the flute geometry and the materials used in the fabrication of this instrument. After modelling these losses in the frequency domain, they will be represented using a system approach where the fractional order part is separated from the system’s transfer function. Thus, this representation allows us to study, in a precise way, the influence of the fractional order behaviour on the overall system. Effectively, the fractional behavior only appears much below the 20 Hz audible frequencies, but it explains the influence of this order on the frequency response over the range [20–20,000] Hz. Some simulations will be proposed to show the effects of the fractional order on the system response.


Introduction
Fractional order derivative is an ancient mathematical representation that appeared in the late 16th century, and its actual implementation started fifty years ago. However, the numerical simulation of such systems is relatively new, as special tools need to be integrated in the simulators. Thus, the first efficient and inexpensive simulations introduced for simple models (conservative plane waves) were based on signal processing tools: the so-called digital waveguide formalism [1][2][3][4] and more specifically, a factored form introduced by Kelly-Lochbaum [5]. The initial idea rested on the factoring of the alembertian of the equation of plane waves into two transport operators that each govern decoupled "round trip" progressive waves, from which we could derive a form in efficient delay system for the simulation [6,7].
Furthermore, 3D and 2D models with realistic boundary conditions are far too complex to be considered, for example in real-time sound synthesis. They can be effectively reduced to a 1D wave equation including a term that models the flare of the tube profile. It is about the equation of the pavilions, which is also called model of Webster. A more elaborative version of this conservative model includes the effect of visco-thermal losses due to the boundary layers in the vicinity of the walls. This dissipative model, known as the Webster-Lokshin 1D [8][9][10], includes a term that involves a fractional derivation in time of order 3/2 [11][12][13][14][15]. This operator plays a crucial role from a perceptual point of view on sound realism [16].
Thus, the work presented in this paper is a part of a larger project. The objective is not to study the system from an acoustic musical point of view for which there are indeed many works, but from a control point of view, in particular within the framework of the dynamics of complex systems during the study of the coupling between the nonlinear exciter and the resonator (which is part of the continuity of this paper and which will be the subject of a future publication).
Concerning this work, it consists of modelling and simulating the viscous losses within the resonator of the flute musical instrument. This work is not unique, as it is a continuity of three previous publications that showed the following: • The system consisting of the musician-flute that was implemented and modeled. It consisted of an air compressor, a servo-valve, and an artificial mouth mounted to mimic the musician's lungs and mouth. A control system was also developed to regulate the pressure and the flow delivered to the artificial mouth. Added to that, the flute exciter was directly coupled with the artificial mouth and some pressure and temperature sensors were placed within the resonator [17,18]; • The knowledge model was developed to represent the transfer between the pressure source at the input of the tube at x = 0 and the flow at any point x of the tube of length L and of constant radius r. Partial differential equations aiming to causally decompose the global model into sub-models, and thus to facilitate analysis in the frequency domain, were used in modelling [19].
So, in more detail, the main target of this article is to analyze the effect of viscous thermal fractional order element on the sound delivered at the output of the flute. The novelty of the work resides in the numerical synthesis of the viscous thermal losses as well as in its simulation using the hardware-in-the-loop technique. This numerical simulation is important, as it allows the sweeping of the fractional order viscous thermal variable, which is mainly linked to the geometry of the flute's resonator as well as to the materials constituting it.
This article will be divided as follows. In Section 2, the modelling of the test bench will be presented. System approach representation as well as the analysis of all the blocks will be introduced in Section 3. Section 4 presents the rationalization technique of a fractional order system in order to represent it using a series of real poles and zeros of an integer order. Finally, Section 5 summarizes this work and proposes some future tasks.

Schematisation, Configuration, and Setting in Equation
Let us consider an acoustic tube of length L and a constant radius r subjected to an acoustic flow (also called volume flow) Qv(t) with Figure 1). When an acoustic wave propagates in the air, this sets the particles of the fluid in motion, which vibrate at a speed v(t) around their equilibrium position. Then, the acous- tic flow Qv(t) measures the flow [in m 3 /s] of this speed through a surface and presents it as a scalar quantity [20][21][22][23][24]. The acoustic impedance Zac (also called specific acoustic impedance, because it is an intensive quantity) of a medium is defined in steady state by the ratio between the acoustic pressure [in Pa] and the speed [in m/s] of the associated particle. When the medium is air, Zac is equal to the product between the density of air, ρa, and the speed of sound in air, ca; thus, Zac = ρa ca. These two parameters depend also on the air temperature Ta. For more illustration, Table 1 gives the values of the speed of the sound ca of the density ρa and of the characteristic acoustic impedance Zac as a function of the temperature Ta of the air. The model used in this work is based on Webster-Lokshin [9]. It is a model with a mono-spatial dependence that characterizes the linear propagation of acoustic waves in tubes with axial symmetry. This model takes also into account viscous-thermal losses at the wall boundaries with the assumption of wide tubes [6]. Thus, in an axisymmetric tube of constant section S = π r 2 , the acoustic pressure P(x,t,L) and the acoustic flow Qv(x,t,L) are governed by the equation of the pavilions, which is also called Webster-Lokshin, and Euler equation, leading to system (1): where ɛ is a parameter associated with visco-thermal losses. More precisely, ɛ is given by the relation: where lv and lh represent the characteristic lengths of viscous (lv = 4 × 10 −8 m) and thermal (lh = 6 × 10 −8 m) effects, γ being the ratio of a specific heat [26].
The phenomenon of visco-thermal losses is a dissipative effect at the wall of the tube, which is due to the viscosity of the air and to the thermal conduction [6,27]. For the case of wind musical instruments resonators, the assumption of wide tubes is used. This hypothesis is expressed by the following relation: where λ = ca/f represents the wavelength (in m) and f the frequency (in Hz). Thus, for a speed of the sound ca and a frequency fmin corresponding to the lower frequency limit of the model to be studied, it is possible to determine the minimum value of the radius rmin of the acoustic tube below which the model is not valid.

Resolution in the Symbolic Domain
Under the assumption of zero initial conditions (I.C = 0), the Laplace transformation applied to the system (1) leads to: , s being the Laplace variable and TL being its transformation.
Solving the Webster-Lokshin equation [9] gives the solution ( ) , , P x s L in the general form: where A(s) and B(s) are rational functions of s that depend on the boundary conditions, and Г(jω) = j k(ω), k(ω) is a standard complex wave number. Г(s) is given in the Laplace domain by the following relation [28]: The expression of the solutions of ( ) is deduced in two steps: • Using the Euler equation in the Laplace domain (second equation of the system (4)), that is, • Introducing the general solution of ( ) in relation (7), that is, Finally, the solution of ( ) Taking into account the boundary conditions makes it possible to determine the two unknowns A(s) and B(s), and finally the impedance ( ) ( ) ( ) = of the finite medium of length L.
From the perspective of a system approach, the function Г(s) defined in relation (6) is rewritten as follows: (10) or again, in canonical form,  (11) where ωr,m is a transitional frequency (in rad/s). This expression is very representative as it allows, from a system approach point of view, the extension of the fractional model to make it possible to easily vary, in numerical simulation, the fractional order m, which is the image of visco-thermal losses, while from an experimental point of view, it would be necessary to manufacture and test a large number of resonators with different dimensions, roughness, and materials.
Note that in the theoretical case where the system is conservative, that is to say ɛ = 0, the function Г(s) (relation 10) is reduced to Г(s) = s/ca. By replacing Г(s) of relation (10) in can be expressed as follows:  (12) or again, by introducing the characteristic acoustic impedance Zac = ρa ca and the transi- becomes: at any point x of the acoustic tube of length L makes it possible to deduce the pressure ( ) [29]. To conclude this paragraph concerning the resolution in the symbolic domain, the study of asymptotic behaviors of Z(x,s,L), that is ω ω (14) and , lim (15) highlights that Z(x, s, L) tends towards a behavior of the type:

-
Fractional derivative of order m + 1, i.e., 1.5 with m = 0.5, when s tends to zero; -Proportional, whose gain value is fixed by Zac/S, when s tends to infinity.

Frequency Response
In a stationary harmonic regime, the frequency response Z(x,jω,L) is given by where the transitional frequencies ωr,m and ωL,x have the following expressions: Thus, ωr,m decreases when the radius r increases, and on the contrary, ωL,x increases when the position x moves away from the origin and approaches the end L of the acoustic tube.

System Approach
From a causal point of view, the input of the resonator at x = 0 is defined by the pressure at the output of the nonlinear exciter. This is the reason why the system approach developed in this paragraph considers the admittance Y(x,s,L) = Z −1 (x,s,L) and not the impedance Z(x,s,L). More specifically, it is the input admittance at x = 0, which is denoted Yin(s,L) = Y(0,s,L). Note that this consideration of the admittance Yin(jω,L) leads to an integrative behavior for frequencies lower than ωL,x (derivative for Z(x,s,L)), thus respecting integral causality, which is a fundamental notion in a system approach.
In addition, the admittance Y(x,s,L) is broken down into a cascade of local transfer functions of which all the parameters, as well as all the input and output variables, have a physical meaning. Then, this decomposition facilitates the frequency analysis of the Webster-Lokshin model, thus reaching a reduced model to be implemented in the simulator.

Decomposition of Admittance Y(x, s, L) into Subsystems
Therefore, the admittance Y(x,s,L) = Z −1 (x,s,L) of an acoustic tube of length L at a point x between 0 and L is defined by the expression: ω ω ω ω ω (18) relation that can be expressed as follows: where, For the following, the concept of acoustic admittance , which is defined between the pressure source At x = 0, for this finite medium, the admittance of input Yin(s,L) therefore has the expression: Always for x = 0, but for a semi-infinite medium ( Figure 2 presents the block diagrams associated with this system approach where the different transfer functions are defined by: (c) Figure 2. Block diagrams associated with the system approach: whenever x is between 0 and L (a), at x = 0 for the finite system L (b), at x = 0 for a semi finite system (c).

Frequency Analysis of the System Approach
In a stationary harmonic regime, the relation (24) becomes: The remaining part of this paragraph is devoted to a detailed analysis of the frequency responses Im(jω), F(x,jω,L), and T(x,jω,L) of each subsystem. An analysis of the whole system response H(x,jω,L) will be also discussed.

Analysis of Im(jω)
The analysis of Im(jω) highlights two behaviors whose transition zone is fixed by the transitional frequency ωr,m; these behaviors are: • For ω << ωr,m, a fractional integrative behavior of order m/2 = 0.25. Indeed, The values of r = 5 mm and L = 0.3 m correspond to the dimensions of the experimental device developed in a first part of the overall project and used to validate a numerical simulator of the artificial mouth + nonlinear exciter + resonator assembly, in addition to a simulator developed using MatLab/Simulink.
Thus, for the area of study considered in this work, the area defined by the range [20-20,000] Hz of the audible frequencies, the transfer Im(s) can be reduced to the unit, which amounts to writing , allowing a reduction in the block diagrams of Figure 2. The direct consequence is that in the case of a semi-infinite medium at x = 0, the fractional integration behavior has no influence on the range of audible frequencies.  Knowing that in the case of a recorder ωr,m << ωL,x, the analysis of F(0,jω,L) again highlights two behaviors whose transition zone is fixed by the transitional frequency ωr,m, that is:

Analysis of T(0,jω,L)
The analysis of T(x,jω,L) highlights three behaviors whose transition zones are fixed by the transitional frequencies ωr,m and ωL,x: • For w << ωL,x, an integrative behavior with two different orders according to the frequency range. Indeed,  and, for ωr,m << ω, a derivative behavior of order 1, that is

Analysis of H(x,jω,L)
The analysis of H(x,jω,L) highlights three behaviors whose transition zones are fixed by the transitional frequencies ωr and ωL,x: • For ω << ωr << ωL,x, an orderly fractional integrative behavior -(1 − m/2) = −0.75, that is • For ωr,m << ω << ωL,x, a derivative behavior of order 1, that is For ωL,x << ω, a behavior composed of an alternation of anti-resonances and resonances, that is, Note that the farther the position x moves away from the origin, the higher the transitional frequency ωL,x pushes the anti-resonance and resonance frequencies toward the high frequencies.
In addition, the position x having no influence on the transitional frequency ωr,m, the fractional integrative behavior of order −0.75 still does not appear in this frequency range.

Study of the Influence of the Fractional Order m
In the Webster-Lokshin model, the fractional order m has the value 0.5. The objective of this paragraph is to analyze the influence of the order m on the behavior of the resonator by considering that m belongs to the interval [0; 1] with a nominal value m0 = 0.5, which is a consideration that facilitates the introduction of the concept of parametric uncertainty (additive or multiplicative) at the fractional order level. Thus, by generalizing the expression of the parameter ε = K0/r associated with visco-thermal losses in the Webster-Lokshin model at ε = 2 m K0/r which for m = 0.5 gives the same expression), the analytical link is naturally established between visco-thermal losses and fractional order.
Thus, the fractional order occurs only in the presence of visco-thermal losses. In the theoretical case of a purely conservative system, the parameter ε is zero, which is equivalent to m = 0, taking into account the relation (44). In this case, the expression of the acoustic transfer H (x, s, L), denoted then H0 (x, s, L), of a finite medium is reduced to: (45) Figure 8 shows the Bode diagrams at x = 0 of H(0,jω,L) for different values of the fractional order over the range  Hz of the audible and achievable frequencies with a recorder (Figure 7). In order to magnify the different curves in Figure 8 for better observation, Figure 9 presents the reduced frequency responses H(0,jω,L)/H0 with the frequency axis on a linear scale over the range  Hz.
The observation of these frequency responses shows that the influence of the order m is essentially located: • For gain diagrams, at the peaks of resonances and anti-resonances; quantifiable effects using quality factors Qzi for anti-resonances and Qpi for resonances clearly illustrate the phenomenon of dissipation associated with visco-thermal losses.
• For phase diagrams, at the crossing points at 0° with a local slope, which is all the more important as the order is small, and the slope becomes infinite for m = 0 (purely conservative case).

From the Simplified Fractional Model to its Rational Forms
For the area of study defined by the range [20-20,000] Hz of the audible frequencies, the analysis presented in the previous paragraph shows that the frequency response Im(jω) is: which can be reduced to the unit. This is the reason why for this field of study, the frequency response H(x,jω,L) defined, as a reminder, by is simplified and noted as ( ) ω , that is: with, always as a reminder, (49) Figure 10 shows the block diagrams associated with this simplification in the field of study. In general, the temporal simulation of fractional models often requires the use of rational models [30]. Thus, the fractional form of defined by the relation (49) can be put in a rational form of N cells in cascade, noted  (50) to the parallel form (51) by decomposing into simple elements. Note that the parallel rational form facilitates the return to the time domain by inverse Laplace transform, and that it is often associated with a decomposition in a modal space [31].
From a theoretical point of view, the ωzi corresponds to the roots of the numerator of T(x,jω,L), that is: From a practical point of view, finding these roots by analytical resolution is complex, if not impossible. On the other hand, the search by numerical resolution does not pose any particular problem. For example, it is possible to use the fact that the alternation of ωzi and ωpi appears clearly on the phase from when passing at 0° from − 90° to + 90° for ωzi and from + 90° to − 90° for ωpi (see example illustration below).
In the context of the work of this thesis, the rational forms of N cells in cascade and in parallel are considered as behavior models whose numerical values of the parameters are obtained using an optimal approach aiming to minimize the difference between the This digital procedure is available in the Frequency Domain System Identification (FDSI) module of the CRONE Toolbox [32].
The procedure consists, in a first step, in generating the target frequency response ( ) ω of the fractional form, which appears in red (Data) in Figure 11a. Then, cells are added one after the other by clicking on the "Add New Cell" command in the "Cells actions" menu (in green at the top right), then by positioning the mouse cursor on the phase diagram at a point considered where the cutting phase of ( ) axis 0°, and going from the lowest values (left) to the most important (right). Note that in this graphical interface, the term "cell" corresponds to a polynomial (numerator or denominator). Thus, with each addition, a column in blue appears in the "Transfer function" menu (in purple at the bottom) for the numerator and in a column in yellow appears for the denominator. The first line "Cell Frequency" gives the value in rad/s of ωzi or ωpi, the second "Cell Order" gives the highest order of the polynomial (here +2 for the numerator and −2 for the denominator), the third "Local Order" is equal to +1 (for the numerator) and −1 (for the denominator) insofar as these are explicit forms [27], and finally, the fourth "Damping Factor" gives the value of ζzi or ζpi. All these values in the columns can be modified at will by clicking in the corresponding box. Thus, in the case of the resonator, all the values of ζzi or ζpi are initialized to 0.01. Therefore, this first stage of the procedure makes it possible to fix the structure of the behavior model, as well as the initial values of its parameters. The second step is an optimization step launched using the "Tools" menu (in orange at the bottom right). For the example of illustration, the result appears in (Figure 11b) with in particular the optimal values of ωzi, ωpi, ζzi, and ζpi.
Thus, the blue curve corresponds to the frequency response

Remark 1.
In linear systems dynamics, in the general case of a polynomial of order 2 having one pair of conjugate complex roots, the damping factor ζ and the resonance factor Q associated with this pair are linked by a relation of the form: In instrument acoustics [28] and in particular in the specific case of resonators of wind instruments, the damping factors are very small in front of the unit. This is the reason why the relation (54) is reduced to: Note that in instrument acoustics, the term quality factor is used in place of the resonance factor. Thus, in many works, taking the visco-thermal losses into account is made directly using the parallel rational form defined by the relation (54) in which the ζpi are replaced by the corresponding Qpi (relation (55)) [26] without going through fractional models. Table 2 summarizes the final numerical values of the parameters ωzi, ζzi, Qzi, ωpi, ζpi, and Qpi of the N = 4 cells of the cascade form to which we must not forget the cell number zero, namely the integrator A0/s. As for Table 3, it gives the numerical values of the parameters Ai, Bi, ωpi, ζpi, and Qpi of N = 4 cells of the parallel form (relation (55)) to which is added cell A0/s.

Conclusions and Future Works
The structure and progression of this article are organized in a didactic way so that readers with no idea about visco-thermal losses in wind instruments can "absorb" the dynamic behavior of an acoustic tube of constant radius. From the two partial differential equations that define the Webster-Lokshin model, a classical resolution in the operational domain leads to the analytical expression of the acoustic impedance and admittance of the function tube of position x, its length L, and its radius r. Moreover, a system vision is proposed aiming to causally decompose the global model into sub-models, thus facilitating analysis in the frequency domain. One of the conclusions of this frequency analysis is that the fractional model can be simplified over the range 000] Hz of the audible frequencies. In addition, the introduction of an uncertainty at the level of the fractional order (whose value considered as nominal is that of the initial Webster-Lokshin model, namely m0 = 0.5) allows us to study the influence of the order m when this varies between 0 (conservative case) and 1.
As is often the case with fractional models, simulation in the time domain requires the establishment of rational forms. Thus, two rational forms composed of an integrator and N second-order cells, one in cascade and the other in parallel, were introduced. Then, the parameters of the cascade form are determined using the Frequency Domain System Identification (FDSI) module of the CRONE Toolbox. As for the parameters of the parallel form, they are obtained by a decomposition into simple elements of the cascade form.
More generally in the fractional model, this study of visco-thermal losses within the resonator of a wind instrument leads to a finding similar to that already made in other fields. Indeed, the main interest of the fractional form resides in the parametric parsimony, that is to say, the capacity that the integro-differential operator of non-integer order has to model with a minimum of parameters the greatest number dynamic phenomena. Thus, the study of parametric sensitivity, in particular in the frequency domain, is simpler.
As a future work, a comparison between the simulated values and the exact outputs can be computed in order to compare the approximation effects from a practical point of view. Thus, the model uncertainties will be analyzed in more detail when comparing the real and the simulated systems. Added to that, building resonators with the same fractional order as proposed in this article will be a major challenge, as going from simulated systems to implemented ones will be an innovation in the musical instruments field. So, in more detail, an extension of the fractional model to take into account the visco-elastic losses is proposed, thus making it possible to vary the fractional order m from 0 (conservative system) to 1 (very dissipative system), and not to consider only m = 0.5 as is currently the case in the literature. This domain [0; 1] belonging to the order m makes it possible to better account for the influence of geometry (radius r and length L), the roughness, and the nature of the material of the resonator.