© The Author(s), Published by EDP Sciences, 2023

1 Introduction

Sound synthesis by modal decomposition is based on an input impedance measurement, which captures the passive acoustical response of a real instrument. This method thus has the advantage of capturing the acoustical subtleties that differentiate two clarinets by means of a single measurement. Other methods, such as delay lines [1], waveguides [2, 3] or spatial discretization [4, 5] would require more effort on the geometrical description of the resonator to render such acoustical subtleties.

Furthermore, sound synthesis by modal decomposition requires little RAM compared to such other methods, which can be beneficial for embedded devices. Indeed, in modal decomposition synthesis, the temporal integration scheme requires to keep only a few previous iterations (precisely two for [1]) to compute a new one. For the delay line synthesis, it is necessary to keep in memory all the temporal iterations during a round trip, i.e. during 2LF_s/c₀ iterations for an academic cylindrical resonator without lateral holes. Although this is not a problem for most embedded processors, fine modeling of viscothermal losses in the waveguide formalism also requires additional computations involving fractional derivatives [6]. Finally, the ODE formalism is particularly well suited to the bifurcation analysis, which is not the case of the PDE [6] or of the DDE formalisms [7].

The input impedance is a measure of the linear frequency response of the resonator, for a low-amplitude excitation. However, in the case of woodwind instrument playing, the measured acoustic pressure and velocity in the resonator can be very high. For instance, [8] measure a pressure in the mouthpiece of 4 kPa, i.e. 163 dB for a monochromatic wave. According to [9] of Chap. 8.4.5, from an acoustic speed of roughly 1 m/s, jet separation phenomena appear. At the level of a geometrical discontinuity, such as the open end of a pipe [10] or a lateral orifice [11], vortices are observed. A part of the kinetic energy of the jet is absorbed by these vortices and is dissipated as heat by friction. These losses are accounted for through a nonlinear relationship derived from Bernoulli’s law, as demonstrated by the implementation in waveguide modeling or delay lines of Refs. [7, 12–14].

Following [11], the works from Atig et al. [12, 15] and Dalmont et al. [13, 16] model localized nonlinear losses through a resistive impedance Z_t. This impedance is in series with the radiation impedance Z_R, and depends on the amplitude of the acoustic velocity v₀ at the location of the discontinuity:

$$Z_t=\frac{4c_d}{3\pi }\frac{\left|v_0\right|}{c_0}{Z}_c,$$(1)where Z_c = ρ₀c₀/S is the characteristic impedance of the medium for plane waves, and c_d ∈ [0; 3] is a parameter depending on the geometry of the termination. The larger the radius of curvature at the output, the smaller this coefficient is. The coefficient c_d is difficult to predict theoretically, and easier to determine experimentally [12]. As shown by Atig et al. [12], nonlinear losses have a significant influence on the playing range of a clarinet, hence models of sound production in woodwinds should include such effects.

Models of woodwind instruments taking into account localized nonlinear losses at the end of the resonator have never been developped in the framework of modal synthesis. This is the main contribution of this work.

First, a method proposed by Diab et al. [17] to account for nonlinear losses by modal decomposition is presented. It is then applied to a cylindrical tube with nonlinear losses located at the open end. From this resonator, a self-oscillating clarinet-like system is defined and simulated by time integration. These simulations are finally compared to the experimental results published by Atig et al. [12] and Dalmont and Frappé [13].

2 Method of modal decomposition accounting for localized nonlinear losses

The definition of a “nonlinear impedance” (i.e. adjusting the impedance so that the link between acoustic flow and pressure stays valid in nonlinear conditions) opens the way to consider localized nonlinear losses by modal decomposition. In Diab et al. [17], the evolution of the surface impedance of perforated plates is computed by temporal simulation for a broadband excitation, with respect to the RMS amplitude of the acoustic velocity at the level of the hole, noted v_RMS. From the decomposition of the impedance in the linear domain into a sum of N modes, two methods are studied for increasing values of v_RMS: the interpolation of the impedance, on one hand, and the regression of the poles and complex residues (s_n, C_n) on the other hand. This second method will be considered hereafter.

Nonlinear effects are taken into account by allowing the poles and residues to vary with respect to v_RMS. In the Laplace domain, the modal decomposition of the input impedance is written:

$$Z(s,v_{\mathrm{RMS}})=Z_c\sum _{n=1}^N\mathrm{}\frac{C_n\left(v_{\mathrm{RMS}}\right)}{s-s_n\left(v_{\mathrm{RMS}}\right)}+\frac{C_n^{\mathrm{*}}\left(v_{\mathrm{RMS}}\right)}{s-s_n^{\mathrm{*}}\left(v_{\mathrm{RMS}}\right)},$$(2)where s is the Laplace variable and •* denotes the complex conjugation operation. The relationship between poles and residues and v_RMS is assumed to be a rational fraction [17]. In the present work, the approximation is limited to a polynomial of degree N_p:

$$C_n\left(v_{\mathrm{RMS}}\right)=C_n^{\left(0\right)}+\sum _{k=1}^{N_p}\mathrm{}{C}_n^{\left(k\right)}{\left(v_{\mathrm{RMS}}\right)}^k,$$(3)

$$s_n\left(v_{\mathrm{RMS}}\right)=s_n^{\left(0\right)}+\sum _{k=1}^{N_p}\mathrm{}{s}_n^{\left(k\right)}{\left(v_{\mathrm{RMS}}\right)}^k,$$(4)where the index •⁽⁰⁾ refers to the linear part of the modal coefficient. In Diab et al. [17], the regression coefficients of the poles and residues are then adjusted on a nonlinear surface impedance model Z_NL (v_RMS), given by Laly et al. [18]. It is worth noting that in the present paper, the poles and residues are fitted on the dimensioned quantity of the RMS velocity. Therefore, the unit of coefficients $C_n^{\left(k\right)}$ and $s_n^{\left(k\right)}$ are expressed in $\mathrm{rad}\cdot {\mathrm{s}}^{-1}\cdot [\mathrm{m}\cdot {\mathrm{s}}^{-1}{]}^{-k}$.

The last step consists in computing v_RMS, given by:

$$v_{\mathrm{RMS}}^2\left(t\right)=\frac{1}{t}{\int }_0^t\mathrm{}{v_0}^2\left(\tau \right)\mathrm{d}\tau ,$$(5)where v₀ (t) = u₀/S is the acoustic velocity at the discontinuity, of cross section S. This definition is rather impractical in the framework of a time-domain simulation. In order to obtain a state formalism, it is replaced by:

$$\frac{\mathrm{\partial }\left(t{v}_{\mathrm{RMS}}^2\right)}{\mathrm{\partial }t}=v_0(t{)}^2,$$(6)which is equivalent. By considering the right-hand side of equation (6) as an input for the ODE solver, $t{v}_{\mathrm{RMS}}^2$ can be computed, hence v_RMS.

3 Application to a cylindrical tube with nonlinear losses at the open end

The impedance regression technique presented by Diab et al. [17] is applied to the case of a closed-open cylinder of radius R and length L, which is a similar case study to Ref. [12]. Viscothermal losses are taken into account in the propagation:

$$\mathrm{\Gamma }\left(s\right)=\frac{s}{c_0}+(1+j)\frac{\eta \sqrt{s}}{R\sqrt{2j\pi }},$$(7)where η = 3 × 10⁻⁵ s^1/2. The boundary condition at x = L is given by the radiation impedance Z_R, which includes in series the nonlinear impedance Z_t defined by equation (1):

$$\begin{array}{ll}{Z}_R& =Z_R^{\left(\mathrm{lin}\right)}+Z_t,\\ \mathrm{where} Z_R^{\left(\mathrm{lin}\right)}& =Z_c\left(\mathrm{jk}\Delta l+\frac{1}{4}(\mathrm{kR}{)}^2\right),\\ Z_t& =Z_c\frac{v_{\mathrm{RMS}}(L,t)}{c_0}{K}_{\mathrm{NL}},\end{array}$$(8)and jk = s/c₀, Δl ≈ 0.6R ([9], Chap. 12.6.1.3), K_NL = 4c_d/(3π). In the following, v_RMS(L, t) will be denoted as v_RMS for reading comfort. The input impedance is defined by:

$$z_{\mathrm{in}}=\tanh \left(\mathrm{\Gamma }L+{\tanh }^{-1}\left(z_R\right)\right),$$(9)introducing the dimensionless notation z_• = Z_•/Z_c. A linear input impedance $z_{\mathrm{in}}^{\left(\mathrm{lin}\right)}$ is also defined, such that

$$z_{\mathrm{in}}^{\left(\mathrm{lin}\right)}=\tanh \left(\mathrm{\Gamma }L+{\tanh }^{-1}\left[\mathrm{jk}\Delta l+\frac{1}{4}(\mathrm{kR}{)}^2\right]\right),$$(10)which is equivalent to z_in for v_RMS = 0.

Figure 1 illustrates the evolution of z_in (computed through equation (9) and equation (10), in which the Laplace variable is substituted by j2πf) as a function of the RMS velocity at the open end. It can be observed that taking into account nonlinear losses has a consequence on the amplitude, but not on the frequency of the resonance peaks. Figure 1b shows that this impact is particularly accentuated on the first resonance peak: its amplitude is 43% lower for v_RMS = 24 m/s than in the case without nonlinear losses. The difference is reduced to 30% for the second peak, and to 16% for the sixth peak. Near the antiresonances, the negative relative discrepancies are high, since the impedance values are all close to zero. The absolute deviation remains very low. From the lack of frequency variation of the resonance peaks with v_RMS, we could anticipate that nonlinear losses at the end of the pipe would rather impact dynamics than intonation.

Figure 1 Input impedance of a cylinder zin (computed through Eq. (9) and Eq. (10)) of length L = 64 cm and radius R = 8 mm for different values of the acoustic RMS velocity at the open end. The value of cd has been set to 13/9, which corresponds to a radius of curvature of 0.3 mm, according to [29]. From top to bottom, from left to right: modulus of the input impedance; relative gap between nonlinear and linear definition of zin; detail view on the first peak; detail view on the fourth peak.

3.1 Modal decomposition of the nonlinear input impedance

In order to take into account nonlinear losses in temporal simulation, the modal coefficients of the input impedance must be determined. In the same way that the nonlinear input impedance depends on the acoustic velocity at the open end, the coefficients C_n and s_n derived from the modal decomposition of z_in also depend on v_RMS. To formulate the relationship between the modal coefficients and v_RMS, a minor adaptation of the method presented in Ref. [7] of Chap. 5.5.3 is used. The definition of z_R used in [9] is substituted by the expression given by equation (8), taking into account nonlinear losses at the end of the tube.

3.1.1 Expression of the poles s_n

Equation (9) can also be written as:

$$\begin{array}{l}{z}_{\mathrm{in}}=\frac{\sinh \left[\mathrm{\Gamma }\left(s\right)L+h(s,v_{\mathrm{RMS}})\right]}{\cosh \left[\mathrm{\Gamma }\left(s\right)L+h(s,v_{\mathrm{RMS}})\right]},\\ \mathrm{where} h(s,v_{\mathrm{RMS}})={\tanh }^{-1}\left(z_R\right).\end{array}$$(11)

The poles $s_n\left(v_{\mathrm{RMS}}\right)$ are the solutions of cosh[ΓL + h] = 0, i.e.:

$$\mathrm{\Gamma }\left(s_n\right)L+h\left(s_n,v_{\mathrm{RMS}}\right)-j\frac{\left(2n-1\right)\pi }{2}=0.$$(12)

The solutions of equation (12) give s_n, for different input values of v_RMS. They are plotted on Figure 2a. The consideration of nonlinear losses at the end of a cylindrical pipe has almost no influence on $\mathfrak{I}\left(s_n\right)$, which is the resonance angular frequency of peak n. However, a larger $v_{\mathrm{RMS}}$ increases $\left|\mathfrak{R}\right(s_n\left)\right|$. This shift remains almost the same regardless of the index of the pole. However, in terms of relative deviation, the ratio $\left|\mathfrak{R}\right(s_n-s_n^{\left(\mathrm{lin}\right)})/\mathfrak{I}(s_n\left)\right|$ is decreasing as n increases. This is signaled by the amplitude of peaks of higher index appearing less affected by nonlinear losses, as shown in Figure 1b.

Figure 2 Graphical representation of the N = 8 first poles and residues. The values of vRMS are linearly chosen between 0 m/s and 24 m/s, according to [29]. As in Figure 1, cd = 13/9. (a) Real and imaginary parts of poles sn (vRMS), calculated by solving equation (12). Detailed view on s1. (b) Real and imaginary parts of residues Cn (vRMS), calculated with equation (14). Detailed views on C1 and C4.

3.1.2 Expression of the residues C_n

It remains to calculate residue C_n (v_RMS). In the vicinity of a pole s_n, the denominator D of z_in can be written through a first-order series expansion of cosh:

$$\begin{array}{ll}D& =(s-s_n)D\mathrm{\prime }(s_n,v_{\mathrm{RMS}})\\ \mathrm{where} D\mathrm{\prime }& =(\mathrm{\Gamma }\mathrm{\text{'}}L+h\mathrm{\prime })\sinh \left(\mathrm{\Gamma }L+h\right),\end{array}$$(13)introducing notation •′ = ∂•/∂s. According to Chaigne and Kergomard ([7], Chap. 5.5.3), after application of the residues theorem, the modal decomposition of the input impedance for the nth peak can be written as

$$\begin{array}{ll}{Z}_{\mathrm{in},n}& =\frac{P_n\left(s\right)}{U\left(s\right)}=Z_c\frac{C_n}{s-s_n},\mathrm{where}\\ C_n\left(v_{\mathrm{RMS}}\right)& =\frac{1}{\mathrm{\Gamma }\mathrm{\text{'}}\left(s_n\right)L+h\mathrm{\prime }(s_n,v_{\mathrm{RMS}})}.\end{array}$$(14)

The evolution of coefficients C_n with respect to v_RMS is plotted on Figure 2b. Real and imaginary parts of C_n slightly decrease when v_RMS increases. Following Diab et al. [17], it remains to fit the modal coefficients with respect to v_RMS, as in the example of equation (3).

3.1.3 Polynomial fitting of the modal coefficients

Figure 3a shows the mean relative fitting error for polynomials of different degrees, for the same data as in Figure 2. For linear regression, the mean relative error stays around 10⁻⁷ for each mode n, both for C_n and s_n. Further simulations reveal that choosing an excessive values of c_d tends to increase this error. For a tube with sharp edges at the open end, Atig et al. [12] estimated a maximum value of c_d = 2.8. For c_d = 5, for instance, the mean relative error of linear regression is around 10⁻⁵. Although this error has increased, it remains very low.

Figure 3 Error due to the regression of modal coefficients Cn and sn (a), and its consequence on the discrepancy between modal decomposition and definition (9) on zin (b). (a) Mean relative error ε between fitted modal coefficient $\tilde{X}=\{{\tilde{s}}_n,{\tilde{C}}_n\}$ and its actual value X = {sn, Cn}. Regression is based on Nv = 25 values of vRMS. (b) L2-norm of the error over zin for every frequencies, between definition $z_{\mathrm{in}}^{\left(\mathrm{def}\right)}$ given by equation (9) and modal decomposition approximation $z_{\mathrm{in}}^{\left(\mathrm{modal}\right)}$ defined by equation (2). This error is plotted with respect to the degree of the fitting polynomial of Cn and sn.

Moreover, Figure 3b shows that choosing an polynomial regression degree larger than 1 has no consequence on the error over z_in between modal decomposition with fitted coefficients and the definition given by equation (9). This error is therefore inherent to the modal decomposition approximation (Eq. (2)) of the input impedance (Eq. (9)).

According to these previous results, in the rest of this article, the modal coefficients will both be fitted by a polynomial of order 1.

4 Application to the time-integration simulation of a clarinet-like system

The cylindrical tube considered in Section 3 is now treated as the resonator of a clarinet-like instrument. The physical model governing the self-sustained oscillations of the instrument is presented in Section 4.1. It is then simulated by time integration, imposing a linear increase of the blowing pressure (crescendo) followed by a linear decrease (diminuendo). Details and results from the simulations are given in Section 4.2.

4.1 Physical model for the self-sustained oscillations

The self-sustained system consists of three main parts: reed dynamics, reed channel and resonator. This last subsystem is enriched compared to the classical modal decomposition formalism [19–21] by the addition of nonlinear losses at the end of the pipe.

4.1.1 Reed dynamics

The dynamic behavior of the reed is modeled for its dimensionless displacement x(t) by the following equation:

$$\frac{1}{{\omega }_r^2}\ddot{x}\left(t\right)+\frac{q_r}{{\omega }_r}\dot{x}\left(t\right)+x\left(t\right)=p\left(t\right)-\gamma \left(t\right)+F_c(x,\dot{x}),$$(15)where ω_r = 2π × 2200 rad/s is the reed resonance angular frequency [22] and q_r = 0.4 is the reed damping [23]. γ(t) = P_m(t)/p_M is the dimensionless blowing parameter and p(t) = p^dim (t)/p_M is the dimensionless pressure at the input of the resonator. The beating pressure p_M is the value of the blowing pressure P_m for which the reed closes the reed channel, in quasistatic regime. $F_c(x,\dot{x})$ is the contact force function. In the following, a model of “ghost reed” [24] is considered, i.e. F_c = 0.

4.1.2 Reed channel

The characteristic of the flow u(t) through the reed channel is described by the following equation, according to [25]:

$$u\left(t\right)=-\lambda \dot{x}\left(t\right)+\zeta {\left[x\left(t\right)+1\right]}^{+}\mathrm{sgn}\left[\gamma \right(t)-p(t\left)\right]\sqrt{\left|\gamma \left(t\right)-p\left(t\right)\right|},$$(16)where λ = 5.5 × 10⁻³/c₀ is the reed flow parameter [26], ζ is the embouchure parameter, and γ(t) is the dimensionless blowing pressure. The operator [•]⁺ = (• + |•|)/2 refers to the positive part function. For the simulations, the absolute value functions are regularized to smooth the irregularities of the quasi-static flow characteristic, without altering its behavior, according to [27]. Thus, $|•|\mapsto \sqrt{{•}^2+\eta }$, with η = 0.001.

4.1.3 Resonator

The resonator is described in the frequency domain by its input impedance z_in, under its modal decomposition form given by equation (2). The modal poles and residues s_n and C_n are chosen to be linearly dependent on v_RMS, i.e.:

$$C_n\left(v_{\mathrm{RMS}}\right)=C_n^{\left(0\right)}+C_n^{\left(1\right)}{v}_{\mathrm{RMS}},$$

$$s_n\left(v_{\mathrm{RMS}}\right)=s_n^{\left(0\right)}+s_n^{\left(1\right)}{v}_{\mathrm{RMS}}.$$

The modal decomposition of z_in allows to write the relation between the dimensionless pressure p and flow u in the temporal domain:

$${\dot{p}}_n\left(t\right)-s_n\left(v_{\mathrm{RMS}}\right)p_n\left(t\right)=C_n\left(v_{\mathrm{RMS}}\right)u\left(t\right),$$

$$\mathrm{where}\hspace{1em}p\left(t\right)=2\sum _n^N\mathrm{}\mathfrak{R}\left(p_n\right(t\left)\right).$$(17)

The values of s_n and C_n must be updated at each time step by computing v_RMS (L, t). The method for calculating the mean velocity at nonlinearity is detailed in the following section.

4.1.4 Computation of v_RMS at each time step

To compute v_RMS (L, t), the pressure at the termination L must be calculated. First, a linear problem is considered. The value of s_n is therefore chosen for v_RMS = 0 m/s. The modal components of the pressure are related to the total pressure at L through the following equation:

$$\begin{array}{ll}& p(L,t)=2\mathfrak{R}\left(\sum _n^N\mathrm{}{p}_n\left(t\right){\varphi }_n\left(L\right)\right),\\ & \mathrm{where}\hspace{1em}{\varphi }_n\left(\xi \right)=\cosh \left(\mathrm{\Gamma }\right(s_n^{\left(0\right)}\left)\xi \right),\end{array}$$(18)and ξ ∈ [0, L] is the distance along the longitudinal axis of the resonator. For reading comfort, the expression p(t) will be reserved to denote p(0, t).

The dimensionless acoustic velocity at the open end is then calculated from its definition which is related to the pressure field through the dimensionless Euler equation:

$$\frac{\mathrm{\partial }p}{\mathrm{\partial }\xi }(\xi ,t)=-\frac{1}{c_0}\frac{\mathrm{\partial }v}{\mathrm{\partial }t}(\xi ,t).$$(19)

The dimensionless acoustic velocity v(L, t) is thus obtained by numerical integration of the following equation:

$$\dot{v}(L,t)=-2c_0\mathfrak{R}\left(\sum _n^N\mathrm{}{p}_n\left(t\right)\frac{\mathrm{d}{\varphi }_n}{\mathrm{d}\xi }\left(L\right)\right),$$

$$\mathrm{where}\hspace{1em}\frac{\mathrm{d}{\varphi }_n}{\mathrm{d}\xi }\left(\xi \right)=\mathrm{\Gamma }\left(s_n^{\left(0\right)}\right)\sinh \left[\mathrm{\Gamma }\left(s_n^{\left(0\right)}\right)\xi \right].$$(20)

The RMS velocity is finally obtained by double time integration of equation (20), according to equation (6). Since the equations of the system characterizing the self-sustained oscillations are dimensionless, it is necessary to resize v_RMS, because the modal coefficients $C_n^{\left(k\right)}$ and $s_n^{\left(k\right)}$ presented in equation (3) have been regressed from the dimensioned velocity. To do so, $v_{\mathrm{RMS}}^{\left(\mathrm{adim}\right)}$ is multiplied by p_M/(ρ₀c₀). In the following simulations, a value of p_M = 8.5 kPa is set, according to experimental results from [12].

Regarding the calculation of v(L, t), it is not numerically guaranteed that the mean value of $\dot{v}(L,t)$ remains zero. Consequently, the integral $v(L,t)={\int }_0^t\mathrm{}\dot{v}(L,t)\mathrm{d}t$ may diverge. To avoid these divergence problems when integrating $\dot{v}(L,t)$ and v²(L, t), a short-memory term τ is added within each integral. Thus, v(L, t) is now defined by the following convolution product:

$$v(L,t)={\int }_0^t\mathrm{}\dot{v}(L,u)e^{-\frac{t-u}{\tau }}\mathrm{d}u,$$(21)which is written, in the Laplace ($\mathcal{L}$) domain, as:

$$\mathcal{L}\left[v(L,t)\right]=\mathcal{L}\left[\dot{v}(L,t)\right]\frac{1}{s+1/\tau }.$$(22)

By going back to the time domain, equation (21) becomes:

$$\dot{v}(L,t)=-2c_0\mathfrak{R}\left(\sum _n^N\mathrm{}{p}_n\left(t\right)\frac{\mathrm{d}{\mathrm{\varphi }}_n}{\mathrm{d}\xi }\left(L\right)\right)-\frac{1}{\tau }v(L,t).$$(23)

The same method is applied for the computation of v_RMS:

$$v_{\mathrm{RMS}}^2=\frac{1}{\tau }{\int }_0^t\mathrm{}{e}^{-\frac{t-u}{\tau }}{v}^2\left(L,u\right)\mathrm{d}u$$(24)

$$\iff \mathcal{L}\left[v_{\mathrm{RMS}}^2\right]=\frac{1}{\tau }\mathcal{L}\left[v^2(L,t)\right]\frac{1}{s+1/\tau },$$(25)

$$\iff \frac{\mathrm{\partial }\left(\tau v_{\mathrm{RMS}}^2\right)}{\mathrm{\partial }t}=v^2(L,t)-v_{\mathrm{RMS}}^2.$$(26)

4.1.5 Complete system of equations

The complete self-sustained system is governed by the equations defined in Sections 4.1.1, 4.1.2, 4.1.3 and 4.1.4. This article shows the resolution of this problem using an ode solver. The Cauchy problem $\dot{\mathbf{Y}}=\mathcal{F}(\mathbf{Y},t)$ writes, according to equation (28), as:

$$\mathcal{F}\left(\mathbf{Y},t\right)=\mathcal{F}\left(\begin{array}{l}x\left(t\right)\\ \dot{x}\left(t\right)\\ v\left(L,t\right)\\ \tau v_{\mathrm{RMS}}^2\\ p_1\left(t\right)\\ ⋮\\ p_N\left(t\right)\end{array}\right)$$(27)

$$=\left(\begin{array}{l}\dot{x}\left(t\right)\\ -q_r{\omega }_r\dot{x}\left(t\right)+{\omega }_r^2\left[p\left(t\right)-\gamma \left(t\right)-x\left(t\right)\right]\\ -2c_0\mathfrak{R}\left(\sum _n^N\mathrm{}{p}_n\left(t\right)\frac{\mathrm{d}{\mathrm{\varphi }}_n}{\mathrm{d}\xi }\left(L\right)\right)-\frac{1}{\tau }v(L,t)\\ v^2(L,t)-v_{\mathrm{RMS}}^2\\ C_1\left(v_{\mathrm{RMS}}\right)u\left(t\right)+s_1\left(v_{\mathrm{RMS}}\right)p_1\left(t\right)\\ ⋮\\ C_N\left(v_{\mathrm{RMS}}\right)u\left(t\right)+s_N\left(v_{\mathrm{RMS}}\right)p_N\left(t\right)\end{array}\right),$$(28)where u(t) is computed by using equation (16), C_k (v_RMS) are computed with equation (14), and s_k(v_RMS) are computed with equation (12). The initial conditions ${\mathbf{Y}}_{\mathbf{0}}$ will be set to 0 in the following simulations.

4.2 Simulation including nonlinear losses at the open end

The model including nonlinear losses at the open end as represented by the system of equation (28) is now simulated by time integration, using the solver ode45 from Matlab. The input data of the problem are based on experimental results from [12].

4.2.1 Simulation parameters

The parameters related to the resonator, the reed and the reed channel are presented in Table 1. Six values of c_d have been taken from [12], which are between c_d = 0 and c_d = 2.8. The value of ζ = 0.28 was calculated using data corresponding to a “loose embouchure” configuration. For the calculation of $v_{\mathrm{RMS}}$, the short-memory term was set to the period of the first impedance peak, i.e. $\tau =2\pi /\mathfrak{I}\left(s_1\right)$. The evolution of γ is first linearly ascending, from γ = 0 to γ = 3 in 8 s, then linearly descending, from γ = 3 to γ = 0 in 8 s. The total simulation time is therefore 16 s. Simulations were performed for N = 4 modes. A convergence study showed that the dynamic behavior of the system did not change for a higher number of modes.

4.2.2 Results

The bifurcation diagram of the input pressure p^(dim) = p · p_M over the blowing pressure P_m = γ · p_M is represented on Figure 4, with the parameters detailed in Section 4.2.1.

Figure 4 Bifurcation diagram of the L2 norm of p(dim) with respect to Pm, for different values of cd, in an ascending blowing pressure configuration (a) and a descending one (b).

During a crescendo (Fig. 4a), the oscillation threshold is located near 5.7 kPa (γ_th ≈ 0.67). This threshold remains the same, whether nonlinear losses are taken into account (c_d ≠ 0) or not (c_d = 0). However, nonlinear losses have a significant influence on the extinction threshold ${\gamma }_{\mathrm{extup}}$, i.e. the blowing pressure from which the reed stops oscillating and is pressed completely against the mouthpiece. This influence is shown in Table 2. When c_d is increased (i.e. nonlinear losses increase), γ_extup diminishes. Similarly to the experiments conducted by Atig et al. [12], nonlinear losses have an important influence on the dynamic playing range of the musician. Figure 4 shows that when nonlinear losses are low, the range of stable oscillation amplitude that can be obtained is larger.

In the diminuendo phase, the inverse threshold is the same (γ_inv ≈ 0.89) for each geometry at the open end of the pipe. Around this threshold, the amplitude of the input pressure decreases slightly as the losses increase. This behavior is also observed in the experimental curves of Ref. [12] of Figures 12a, 12c and 12e.

This comparaison with the experimental results of Atig et al. is limited to a qualitative study. Parameters such as the reed resonance angular frequency ω_r and reed damping q_r have a strong influence on the dynamics of the system, as demonstrated by Silva et al. [28] for the oscillation threshold. The complete recalibration of the model on experimental results is beyond the scope of this article.

4.3 Comparison with experimental results from Dalmont and Frappé (2007)

An attempt to quantitatively validate the model of nonlinear losses is carried out. To do so, the values of the saturation threshold γ_sat measured by Dalmont and Frappé [13] are employed (see Fig. 5). This threshold corresponds to the blowing pressure γ for which p is maximum during a crescendo. The authors measured this threshold for different reed openings, which are transcribed here in different values of ζ. The dimensioning pressure values p_M used here correspond to the beating pressure estimated by the authors for a diminuendo. Furthermore, the measurements are performed on a cylindrical tube of dimensions L = 50 cm and R = 8 mm. The geometry at the opening has been estimated by the authors at c_d = 2.8.

Figure 5 Experimental results from [13] for the saturation threshold (·), compared with three methods: modal decomposition without nonlinear losses (+); modal decomposition including nonlinear losses (×); analytical solution provided by Raman model including nonlinear losses (*).

The experimental data are compared to simulations using the same parameters L and R as in the experiment. Simulations are performed for two cases: the first one does not take into account nonlinear losses (c_d = 0), the other one takes them into account (c_d = 2.8). The simulations are performed for ascending ramps of blowing pressure varying from $\gamma =0$ to $\gamma =3$ in 10 s. The other parameters are the same as those given in Table 1.

The model presented in this paper is also compared to the Raman model including nonlinear losses proposed by Dalmont and Frappé [13]. In this model, viscothermal losses are simplified by a frequency-independent coefficient α which replaces $\mathfrak{R}\left(\mathrm{\Gamma }\right)$. The simplicity of this model allows the authors to obtain an analytical expression of the saturation threshold ${\gamma }_{\mathrm{sat}}$ as a function of the parameter c_d.

In the present work, the value of α was first adjusted at $\mathfrak{R}\left(\mathrm{\Gamma }\right(j{\omega }_1\left)\right)$, where ${\omega }_1$ is the angular frequency of the impedance peak supporting the oscillation. However, this value produced excessively high estimates of ${\gamma }_{\mathrm{sat}}$ compared to experimental results, especially for high values of $\zeta$. To better match the experimental data, $\alpha =2.7\mathfrak{R}\left(\mathrm{\Gamma }\right(j{\omega }_1\left)\right)\approx 0.147$ was chosen. It should be noted that the analytical expression of ${\gamma }_{\mathrm{sat}}$ exhibits a linear dependence in $1/\alpha$, for low values of α. Thus, ${\gamma }_{\mathrm{sat}}$ is extremely sensitive to the value of α in Raman’s model.

Figure 5 illustrates the evolution of the saturation threshold as a function of $\zeta$ in the experimental case, through modal decomposition simulation, and through the analytical expression based on the Raman model. It appears first that the saturation threshold is overestimated when nonlinear losses are not taken into account (+). In comparison, the two models taking into account nonlinear losses produce results much closer to the experiment. This shows the importance of taking this phenomenon into account in wind instrument simulations.

Furthermore, for the analytical solution given by the Raman model (*), the evolution of ${\gamma }_{\mathrm{sat}}$ follows a straight line whose slope depends on $1/\alpha$. Although the results are close to the experimental ones, the high sensitivity of ${\gamma }_{\mathrm{sat}}$ to α as well as the linear evolution of ${\gamma }_{\mathrm{sat}}$ reflect the excessive simplicity of the Raman model to describe the dynamics of a wind instrument.

Finally, the model presented in this paper (×) has a good agreement with the experimental data, except for very low values of the mouthpiece parameter ($\zeta <0.12$). For these low values of $\zeta$, the model overestimates ${\gamma }_{\mathrm{sat}}$ compared to the experimental results. This overestimation may be related to the transient behavior at extinction caused by the temporal evolution of the control parameter $\gamma$. A characterization of the system by continuation could give results independent of the evolution rate of $\gamma \left(t\right)$.

In conclusion, the correspondence between the results from the model presented in this article and the experimental data from [13] highlights its ability to describe the dynamic behavior of a simplified wind instrument at saturation.

5 Conclusion

This article develops the design of a sound synthesis model of a reed instrument by modal decomposition of the input impedance, taking into account viscothermal losses as well as nonlinear losses at the end of the resonator. The input impedance now depends on the RMS acoustic velocity at a geometric discontinuity (here, the open termination). Poles and residues resulting from the modal decomposition are fitted as a function of this velocity. In a physical model of wind instrument, the pressure-flow relation defined by the resonator is then completed by new equations which account for this dependence with the velocity at the end of the pipe.

To evaluate the ability of the model to reproduce a real phenomenon, comparisons with the experimental results of [12] and [13] have been made. In the first case, simulations show a qualitatively similar behavior regarding the evolution of the extinction threshold depending on the geometry at the open end. In the second case, the model gives a good correspondence with the experimental results, in particular compared to a model without nonlinear losses.

This formalism could be promising in the sound synthesis of wind instruments in real time, which employs the modal decomposition formalism in order to keep a minimum of memory. The inclusion of nonlinear losses in a complete clarinet model could more accurately translate dynamic phenomena essential to the experience of the musician, as suggested by Dalmont and Frappé [13]: “any realistic model of the clarinet should include nonlinear losses in the side holes”.

Conflict of interest

The authors declare no conflict of interest.

Funding

This study has been supported by the French ANR LabCom LIAMFI (ANR-16-LCV2-007-01).

Acknowledgments

The authors thank S. Maugeais and J.-P. Dalmont for their valuable comments.

References

All Tables

All Figures

Figure 1 Input impedance of a cylinder zin (computed through Eq. (9) and Eq. (10)) of length L = 64 cm and radius R = 8 mm for different values of the acoustic RMS velocity at the open end. The value of cd has been set to 13/9, which corresponds to a radius of curvature of 0.3 mm, according to [29]. From top to bottom, from left to right: modulus of the input impedance; relative gap between nonlinear and linear definition of zin; detail view on the first peak; detail view on the fourth peak.

In the text

Figure 2 Graphical representation of the N = 8 first poles and residues. The values of vRMS are linearly chosen between 0 m/s and 24 m/s, according to [29]. As in Figure 1, cd = 13/9. (a) Real and imaginary parts of poles sn (vRMS), calculated by solving equation (12). Detailed view on s1. (b) Real and imaginary parts of residues Cn (vRMS), calculated with equation (14). Detailed views on C1 and C4.

In the text

Figure 3 Error due to the regression of modal coefficients Cn and sn (a), and its consequence on the discrepancy between modal decomposition and definition (9) on zin (b). (a) Mean relative error ε between fitted modal coefficient $\tilde{X}=\{{\tilde{s}}_n,{\tilde{C}}_n\}$ and its actual value X = {sn, Cn}. Regression is based on Nv = 25 values of vRMS. (b) L2-norm of the error over zin for every frequencies, between definition $z_{\mathrm{in}}^{\left(\mathrm{def}\right)}$ given by equation (9) and modal decomposition approximation $z_{\mathrm{in}}^{\left(\mathrm{modal}\right)}$ defined by equation (2). This error is plotted with respect to the degree of the fitting polynomial of Cn and sn.

In the text

	Figure 4 Bifurcation diagram of the L2 norm of p(dim) with respect to Pm, for different values of cd, in an ascending blowing pressure configuration (a) and a descending one (b).
In the text

	Figure 5 Experimental results from [13] for the saturation threshold (·), compared with three methods: modal decomposition without nonlinear losses (+); modal decomposition including nonlinear losses (×); analytical solution provided by Raman model including nonlinear losses (*).
In the text