<figure><img src="image_1.png" /></figure>

HAL
open science

# Amplitude-dependent modal coefficients accounting for localized nonlinear losses in a time-domain integration of woodwind model

Nathan Szwarcberg, Tom Colinot, Christophe Vergez, Michaël Jousserand

► To cite this version:

Nathan Szwarcberg, Tom Colinot, Christophe Vergez, Michaël Jousserand. Amplitude-dependent modal coefficients accounting for localized nonlinear losses in a time-domain integration of woodwind model. Acta Acustica, 2023. hal-03974714v3

HAL Id: hal-03974714

https://hal.science/hal-03974714v3

Submitted on 18 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
# Amplitude-dependent modal coefficients accounting for localized nonlinear losses in a time-domain integration of woodwind model

Nathan Szwarcberg<sup>a, b</sup>, Tom Colinot<sup>a, b</sup>, Christophe Vergez<sup>b</sup>, Michaël Jousserand<sup>a</sup>

## Summary

This article develops the design of a sound synthesis model of a woodwind instrument by modal decomposition of the input impedance, taking into account viscothermal losses as well as localized nonlinear losses at the end of the resonator. This formalism has already been applied by Diab et al. (2022) to the study of forced systems. It is now implemented for self-oscillating systems. The employed method extends the definition of the input impedance to the nonlinear domain by adding a dependence on the RMS acoustic velocity at a geometric discontinuity. The poles and residues resulting from the modal decomposition are fitted as a function of this velocity. Thus, the pressure-flow relation defined by the resonator is completed by new equations which account for the dependence with the velocity at the end of the tube. To assess the ability of the model to reproduce a real phenomenon, comparisons with the experimental results of Atig et al. (2004) and Dalmont et al. (2007) were carried out. Simulations show that the model reproduces these experimental results qualitatively and quantitatively.

## 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 [18], waveguides [23, 22] or spatial discretization [5, 17] 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 [18]) 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_0$ 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 [19].

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

<footnote><sup>a</sup> Buffet Crampon, 5 Rue Maurice Berteaux, 78711 Mantes-la-Ville, France</footnote>

<footnote><sup>b</sup> Aix Marseille Univ, CNRS, Centrale Marseille, LMA, Marseille, France</footnote>
measured acoustic pressure and velocity in the resonator can be very high. For instance, [4] measure a pressure in the mouthpiece of 4 kPa, i.e. 163 dB for a monochromatic wave. According to [7, 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 [15] or a lateral orifice [20], 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 [2], [19], [11] and [28].

Following [20], the works from Atig [2, 3] and Dalmont [13, 11] 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_0$ at the location of the discontinuity:

$$Z_t = \frac{4c_d |v_0|}{3\pi c_0} Z_c, \quad (1)$$

where $Z_c = \rho_0 c_0 / S$ is the characteristic impedance of the medium for plane waves, and $c_d \in [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 [2]. As shown by [2], 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 developed in the framework of modal synthesis. This is the main contribution of this work.

First, a method proposed by [14] 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 [2] and [11].

## 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 [14], 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_{RMS}) = Z_c \sum_{n=1}^{N} \frac{C_n(v_{RMS})}{s - s_n(v_{RMS})} + \frac{C_n^*(v_{RMS})}{s - s_n^*(v_{RMS})}, \quad (2)$$

where $s$ is the Laplace variable and $\bullet^*$ denotes the complex conjugation operation. The relationship between poles and residues and $v_{RMS}$ is assumed to be a rational fraction [14]. In the present work, the approximation is lim-
ited to a polynomial of degree $N_p$:

$$C_n(v_{RMS}) = C_n^{(0)} + \sum_{k=1}^{N_p} C_n^{(k)} (v_{RMS})^k, \quad (3)$$

$$s_n(v_{RMS}) = s_n^{(0)} + \sum_{k=1}^{N_p} s_n^{(k)} (v_{RMS})^k, \quad (4)$$

where the index •(0) refers to the linear part of the modal coefficient. In [14], the regression coefficients of the poles and residues are then adjusted on a nonlinear surface impedance model $Z_{NL}(v_{RMS})$, given by [21]. 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^{(k)}$ and $s_n^{(k)}$ are expressed in rad·s⁻¹·[m·s⁻¹]⁻⁴.

The last step consists in computing $v_{RMS}$, given by:

$$v_{RMS}^2(t) = \frac{1}{t} \int_0^t v_0^2(\tau) d\tau, \quad (5)$$

where $v_0(t) = u_0/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{\partial(tv_{RMS}^2)}{\partial t} = v_0(t)^2, \quad (6)$$

which is equivalent. By considering the right-hand side of Eq. (6) as an input for the ODE solver, $tv_{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 [14] is applied to the case of a closed-open cylinder of radius $R$ and length $L$, which is a similar case study to [2]. Viscothermal losses are taken into account in the propagation:

$$\Gamma(s) = \frac{s}{c_0} + (1 + j) \frac{\eta \sqrt{s}}{R \sqrt{2j\pi}}, \quad (7)$$

where $\eta = 3 \cdot 10^{-5}$ s¹/². 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 Eq. (1):

$$Z_R = Z_R^{(\text{lin})} + Z_t, \\ \text{where } Z_R^{(\text{lin})} = Z_c \left( jk\Delta l + \frac{1}{4}(kR)^2 \right), \quad (8) \\ Z_t = Z_c \frac{v_{RMS}(L,t)}{c_0} K_{\text{NL}},$$

and $jk = s/c_0$, $\Delta l \approx 0.6R$ [7, Chap. 12.6.1.3], $K_{NL} = 4c_d/(3\pi)$. In the following, $v_{RMS}(L,t)$ will be denoted as $v_{RMS}$ for reading comfort. The input impedance is defined by:

$$z_{in} = \tanh \left( \Gamma L + \tanh^{-1}(z_R) \right), \quad (9)$$

introducing the dimensionless notation $z_{\bullet} = Z_{\bullet}/Z_c$. A linear input impedance $z_{in}^{(\text{lin})}$ is also defined, such that

$$z_{in}^{(\text{lin})} = \tanh(\Gamma L + \tanh^{-1}[jk\Delta l + \frac{1}{4}(kR)^2]), \quad (10)$$

which is equivalent to $z_{in}$ for $v_{RMS} = 0$.

Figure 1 illustrates the evolution of $z_{in}$ (computed through Eq. (9) and Eq. (10), in which the Laplace variable is substituted by $j2\pi 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 1 b 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><img src="image_2.png" /><figcaption>Figure 1: Input impedance of a cylinder $z_{in}$ (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 $c_d$ has been set to 13/9, which corresponds to a radius of curvature of 0.3 mm, according to [1]. From top to bottom, from left to right : modulus of the input impedance; relative gap between nonlinear and linear definition of $z_{in}$; detail view on the first peak; detail view on the fourth peak.</figcaption></figure>

## 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 [7, Chap. 5.5.3] is used. The definition of $z_R$ used in [7] is substituted by the expression given by Eq. (8), taking into account nonlinear losses at the end of the tube.

### 3.1.1 Expression of the poles $s_n$

Eq. (9) can also be written as:

$$z_{in} = \frac{\sinh[\Gamma(s)L + h(s, v_{RMS})]}{\cosh[\Gamma(s)L + h(s, v_{RMS})]}, \quad (11)$$

where $h(s, v_{RMS}) = \tanh^{-1}(z_R)$.

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

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

The solutions of Eq. (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 $\Im(s_n)$, which is the resonance angular frequency of peak $n$. However,
a larger $v_{RMS}$ increases $|\Re(s_n)|$. This shift remains almost the same regardless of the index of the pole. However, in terms of relative deviation, the ratio $|\Re(s_n - s_n^{(\text{lin})})/\Im(s_n)|$ 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 1 b.

### 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:

$$D = (s - s_n)D'(s_n, v_{RMS}) \quad (13)$$
where $D' = (\Gamma' L + h') \sinh(\Gamma L + h)$,

introducing notation $\bullet' = \partial \bullet / \partial s$. According to [7, Chap. 5.5.3], after application of the residues theorem, the modal decomposition of the input impedance for the $n$-th peak can be written as

$$Z_{in,n} = \frac{P_n(s)}{U(s)} = Z_c \frac{C_n}{s - s_n}, \quad \text{where} \\ C_n(v_{RMS}) = \frac{1}{\Gamma'(s_n)L + h'(s_n, v_{RMS})}.$$

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 [14], it remains to fit the modal coefficients with respect to $v_{RMS}$, as in the example of Eq. (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^{-7}$ 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, [2] 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^{-5}$. Although this error has increased, it remains very low.

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 Eq. (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 [16, 25, 27] 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
<figure><img src="image_3.png" /><figcaption>(a) Real and imaginary parts of poles $s_n(v_{RMS})$, calculated by solving Eq. (12). Detailed view on $s_1$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(b) Real and imaginary parts of residues $C_n(v_{RMS})$, calculated with Eq. (14). Detailed views on $C_1$ and $C_4$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(a) Mean relative error $\epsilon$ between fitted modal coefficient $\tilde{X} = \{\tilde{s}_n, \tilde{C}_n\}$ and its actual value $X = \{s_n, C_n\}$. Regression is based on $N_v = 25$ values of $v_{RMS}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $L_2$-norm of the error over $z_{in}$ for every frequencies, between definition $z_{in}^{(def)}$ given by Eq. (9) and modal decomposition approximation $z_{in}^{(modal)}$ defined by Eq. (2). This error is plotted with respect to the degree of the fitting polynomial of $C_n$ and $s_n$.</figcaption></figure>

<caption>Figure 3: Error due to the regression of modal coefficients $C_n$ and $s_n$ (Figure 3a), and its consequence on the discrepancy between modal decomposition and definition (9) on $z_{in}$ (Figure 3b).</caption>

following equation :

$$ \frac{1}{\omega_r^2}\ddot{x}(t) + \frac{q_r}{\omega_r}\dot{x}(t) + x(t) = p(t) - \gamma(t) + F_c(x, \dot{x}), \quad (15) $$

where $\omega_r = 2\pi \times 2200$ rad/s is the reed resonance angular frequency [8] and $q_r = 0.4$

is the reed damping [24]. $\gamma(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" [10] 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 [29]:

$$u(t) = -\lambda \dot{x}(t) + \zeta [x(t) + 1]^+ \operatorname{sgn}[\gamma(t) - p(t)] \sqrt{|\gamma(t) - p(t)|}, \quad (16)$$

where $\lambda = 5.5 \cdot 10^{-3}/c_0$ is the reed flow parameter [12], $\zeta$ is the embouchure parameter, and $\gamma(t)$ is the dimensionless blowing pressure. The operator $[\bullet]^+ = (\bullet + |\bullet|)/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 [9]. Thus, $|\bullet| \mapsto \sqrt{\bullet^2 + \eta}$, with $\eta = 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 Eq. (2). The modal poles and residues $s_n$ and $C_n$ are chosen to be linearly dependent on $v_{RMS}$, i.e.:

$$C_n(v_{RMS}) = C_n^{(0)} + C_n^{(1)} v_{RMS}, \\ s_n(v_{RMS}) = s_n^{(0)} + s_n^{(1)} v_{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(t) - s_n(v_{RMS})p_n(t) = C_n(v_{RMS})u(t), \\ \text{where } p(t) = 2 \sum_n \Re(p_n(t)). \quad (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:

$$p(L, t) = 2\Re \left( \sum_n^N p_n(t) \phi_n(L) \right), \quad (18) \\ \text{where } \phi_n(\xi) = \cosh(\Gamma(s_n^{(0)})\xi),$$

and $\xi \in [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{\partial p}{\partial \xi}(\xi, t) = -\frac{1}{c_0} \frac{\partial v}{\partial t}(\xi, t). \quad (19)$$

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

$$\dot{v}(L, t) = -2c_0 \Re \sum_n^N p_n(t) \frac{d\phi_n}{d\xi}(L), \\ \text{where } \frac{d\phi_n}{d\xi}(\xi) = \Gamma(s_n^{(0)}) \sinh[\Gamma(s_n^{(0)})\xi]. \quad (20)$$

The RMS velocity is finally obtained by double time integration of Eq. (20), according to Eq. (6). Since the equations of the system characterizing the self-sustained oscillations are dimensionless, it is necessary to re-size $v_{RMS}$, because the modal coefficients $C_n^{(k)}$ and $s_n^{(k)}$ presented in Eq. (3) have been regressed from the dimensioned velocity. To do so, $v_{RMS}^{(adim)}$ is multiplied by $p_M/(\rho_0 c_0)$. In the following simulations, a value of $p_M = 8.5$ kPa is set, according to experimental results from [2].
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 \dot{v}(L, t)dt$ may diverge. To avoid these divergence problems when integrating $\dot{v}(L, t)$ and $v^2(L, t)$, a short-memory term $\tau$ is added within each integral. Thus, $v(L, t)$ is now defined by the following convolution product:

$$v(L, t) = \int_{0}^{t} \dot{v}(L, u)e^{-\frac{t-u}{\tau}} du, \quad (21)$$

which is written, in the Laplace ($\mathcal{L}$) domain, as:

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

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

$$\dot{v}(L, t) = -2c_0 \Re \left( \sum_{n}^{N} p_n(t) \frac{d\phi_n}{d\xi}(L) \right) - \frac{1}{\tau} v(L, t). \quad (23)$$

The same method is applied for the computation of $v_{RMS}$:

$$v_{RMS}^2 = \frac{1}{\tau} \int_{0}^{t} e^{-\frac{t-u}{\tau}} v^2(L, u) du \quad (24)$$

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

$$\Leftrightarrow \frac{\partial(\tau v_{RMS}^2)}{\partial t} = v^2(L, t) - v_{RMS}^2. \quad (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 Eq. (28), as:

$$\mathcal{F}(\mathbf{Y}, t) = \mathcal{F} \begin{pmatrix} x(t) \\ \dot{x}(t) \\ v(L, t) \\ \tau v_{RMS}^2 \\ p_1(t) \\ \vdots \\ p_N(t) \end{pmatrix} \quad (27)$$

$$= \begin{pmatrix} \dot{x}(t) \\ -q_r \omega_r \dot{x}(t) + \omega_r^2 [p(t) - \gamma(t) - x(t)] \\ -2c_0 \Re \left( \sum_{n}^{N} p_n(t) \frac{d\phi_n}{d\xi}(L) \right) - \frac{1}{\tau} v(L, t) \\ v^2(L, t) - v_{RMS}^2 \\ C_1(v_{RMS})u(t) + s_1(v_{RMS})p_1(t) \\ \vdots \\ C_N(v_{RMS})u(t) + s_N(v_{RMS})p_N(t) \end{pmatrix}, \quad (28)$$

where $u(t)$ is computed by using Eq. (16), $C_k(v_{RMS})$ are computed with Eq. (14), and $s_k(v_{RMS})$ are computed with Eq. (12). The initial conditions $\mathbf{Y}_0$ will be set to $\mathbf{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 Eq. (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 [2].

### 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 [2], which are between $c_d = 0$ and $c_d = 2.8$. The value of $\zeta = 0.28$ was calculated using data corresponding to a "loose embouchure" configuration. For the calculation of $v_{RMS}$, the short-memory term was set to the period of the first impedance peak, i.e. $\tau = 2\pi/\Im(s_1)$. The evolution of $\gamma$ is first linearly ascending, from $\gamma = 0$ to $\gamma = 3$ in 8 s, then linearly descending, from $\gamma = 3$ to $\gamma = 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.
<table><caption>Table 1: Main constants used for the simulations.</caption><thead><tr><th>L (cm)</th><th>R (mm)</th><th>N modes</th><th>ω<sub>r</sub>/2π (Hz)</th></tr></thead><tbody><tr><td>64</td><td>8</td><td>4</td><td>2200</td></tr><tr><th>q<sub>r</sub></th><th>λ (s)</th><th>ζ</th><th>p<sub>M</sub> (kPa)</th></tr><tr><td>0.4</td><td>1.6 · 10<sup>-5</sup></td><td>0.28</td><td>8.5</td></tr></tbody></table>

### 4.2.2 Results

The bifurcation diagram of the input pressure $p^{(dim)} = p \cdot p_M$ over the blowing pressure $P_m = \gamma \cdot p_M$ is represented on Figure 4, with the parameters detailed in section 4.2.1.

<figure><img src="image_7.png" /><figcaption>Figure 4: Bifurcation diagram of the $L_2$ norm of $p^{(dim)}$ with respect to $P_m$, for different values of $c_d$, in an ascending blowing pressure configuration (4a) and a descending one (4b).</figcaption></figure>

<table><caption>Table 2: Values of the extinction threshold $P_{extup} = p_M \cdot \gamma_{extup}$ (kPa), for the simulation of Figure 4a.</caption><tr><td>c<sub>d</sub></td><td>0</td><td>0.15</td><td>0.9</td><td>1.4</td><td>1.7</td><td>2.8</td></tr><tr><td>P<sub>extup</sub></td><td>17.7</td><td>17.3</td><td>15.3</td><td>14.5</td><td>14.2</td><td>13.2</td></tr></table>

During a *crescendo* (Figure 4a), the oscillation threshold is located near 5.7 kPa ($\gamma_{th} \approx 0.67$). This threshold remains the same, whether nonlinear losses are taken into account ($c_d \neq 0$) or not ($c_d = 0$). However, nonlinear losses have a significant influence on the extinction threshold $\gamma_{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), $\gamma_{extup}$ diminishes. Similarly to the experiments conducted by [2], 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 ($\gamma_{inv} \approx 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 [2, Figure 12 a, c and e].

This comparison with the experimental results of Atig et al. is limited to a qualitative study. Parameters such as the reed resonance angular frequency $\omega_r$ and reed damping $q_r$ have a strong influence on the dynamics of the system, as demonstrated by [26] 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 $\gamma_{sat}$ measured by [11] are employed (see Figure 5). This threshold corresponds to the blowing pressure $\gamma$ 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 $\zeta$. 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$.

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 [11]. In this model, viscothermal losses are simplified by a frequency-independent coefficient $\alpha$ which replaces $\Re(\Gamma(j\omega_1))$. The simplicity of this model allows the authors to obtain an analytical expression of the saturation threshold $\gamma_{sat}$ as a function of the parameter $c_d$.

In the present work, the value of $\alpha$ was first adjusted at $\Re(\Gamma(j\omega_1))$, where $\omega_1$ is the angular frequency of the impedance peak supporting the oscillation. However, this value produced excessively high estimates of $\gamma_{sat}$ compared to experimental results, especially for high values of $\zeta$. To better match the experimental data, $\alpha = 2.7\Re(\Gamma(j\omega_1)) \approx 0.147$ was

chosen. It should be noted that the analytical expression of $\gamma_{sat}$ exhibits a linear dependence in $1/\alpha$, for low values of $\alpha$. Thus, $\gamma_{sat}$ is extremely sensitive to the value of $\alpha$ in Raman's model.

<figure><img src="image_8.png" /><figcaption>Figure 5: Experimental results from [11] 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 (*).</figcaption></figure>

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_{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_{sat}$ to $\alpha$ as well as the linear evolution of $\gamma_{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_{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(t)$.

In conclusion, the correspondence between the results from the model presented in this article and the experimental data from [11] 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 [2] and [11] 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 non-

linear 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 [11] : "any realistic model of the clarinet should include nonlinear losses in the side holes".

# Acknowledgments

This study has been supported by the French ANR LabCom LIAMFI (ANR-16-LCV2-007-01). The authors thank S. Maugeais and J.-P. Dalmont for their valuable comments.

# References

[1] Mérouane Atig. *Non-linéarité acoustique localisée à l'extrémité ouverte d'un tube. Mesure, modélisation et application aux instruments à vent*. PhD thesis, Université du Maine, 2004.

[2] Mérouane Atig, Jean-Pierre Dalmont, and Joël Gilbert. Saturation mechanism in clarinet-like instruments, the effect of the localised non-linear losses. *Applied Acoustics*, 65(12):1133-1154, September 2004. Publisher: Elsevier.

[3] Mérouane Atig, Jean-Pierre Dalmont, and Joël Gilbert. Termination impedance of open-ended cylindrical tubes at high sound pressure level. *Comptes Rendus Mécanique*, 332(4):299-304, 2004.

[4] Baptiste Bergeot, André Almeida, Bruno Gazengel, Christophe Vergez, and Diener Ferrand. Response of an artificially blown clarinet to different blowing pressure profiles. *The Journal of the Acoustical Society of America*, 134(4):1199-1204, 2013.
*tical Society of America*, 135(1):479–490, 2014.

[5] Stefan Bilbao. *Numerical sound synthesis: finite difference schemes and simulation in musical acoustics*. John Wiley & Sons, 2009.

[6] Stefan Bilbao and John Chick. Finite difference time domain simulation for the brass instrument bore. *The Journal of the Acoustical Society of America*, 134(5):3860–3871, 2013.

[7] Antoine Chaigne and Jean Kergomard. *Acoustics of musical instruments*. Springer, 2016.

[8] Vasileios Chatziioannou and Maarten Van Walstijn. Estimation of clarinet reed parameters by inverse modelling. *Acta Acustica united with Acustica*, 98(4):629–639, 2012.

[9] Tom Colinot. *Numerical simulation of woodwind dynamics: investigating non-linear sound production behavior in saxophone-like instruments*. PhD thesis, November 2020.

[10] Tom Colinot, Louis Guillot, Christophe Vergez, Philippe Guillemain, Jean-Baptiste Doc, and Bruno Cochelin. Influence of the "ghost reed" simplification on the bifurcation diagram of a saxophone model. *Acta Acustica united with Acustica*, 105(6):1291–1294, 2019.

[11] Jean-Pierre Dalmont and Cyrille Frappé. Oscillation and extinction thresholds of the clarinet: comparison of analytical results and experiments. *The Journal of the Acoustical Society of America*, 122(2):1173–1179, August 2007.

[12] Jean-Pierre Dalmont, Joël Gilbert, and Sébastien Ollivier. Nonlinear characteristics of single-reed instruments: Quasistatic volume flow and reed opening measurements. *The Journal*

Szwarcberg et al., p. 12
*of the Acoustical Society of America*, 114(4):2253–2262, 2003.

[13] Jean-Pierre Dalmont, Cornelis Nederveen, Véronique Dubos, Sébastien Ollivier, Vincent Méserette, and Edwin Sligte. Experimental Determination of the Equivalent Circuit of an Open Side Hole: Linear and Non Linear Behaviour. *Acta Acustica united with Acustica*, 88:567–575, July 2002.

[14] Daher Diab, Didier Dragna, Edouard Salze, and Marie-Annick Galland. Non-linear broadband time-domain admittance boundary condition for duct acoustics. application to perforated plate liners. *Journal of Sound and Vibration*, 528:116892, 2022.

[15] JHM Disselhorst and L Van Wijngaarden. Flow in the exit of open pipes during acoustic resonance. *Journal of Fluid Mechanics*, 99(2):293–319, 1980.

[16] Vincent Dubos, Jean Kergomard, Ali Khettabi, J-P Dalmont, DH Keefe, and CJ Nederveen. Theory of sound propagation in a duct with a branched tube using modal decomposition. *Acta Acustica united with Acustica*, 85(2):153–169, 1999.

[17] Augustin Ernoult, Juliette Chabassier, Samuel Rodriguez, and Augustin Humeau. Full waveform inversion for bore reconstruction of woodwind-like instruments. *Acta Acustica*, 5:47, 2021.

[18] Philippe Guillemain, Jean Kergomard, and Thierry Voinier. Real-time synthesis of clarinet-like instruments using digital impedance models. *The Journal of the Acoustical Society of America*, 118(1):483–494, 2005.

[19] Philippe Guillemain and Jonathan Teroir. Digital synthesis models of clarinet-like instruments including non-linear losses in the resonator. In *of the*
*9th International Conference on Digital Audio Effects*, page 83. Citeseer, 2006.

[20] U Ingård and S Labate. Acoustic circulation effects and the nonlinear impedance of orifices. *The Journal of the Acoustical Society of America*, 22(2):211–218, 1950.

[21] Zacharie Laly, Noureddine Atalla, and Sid-Ali Meslioui. Acoustical modeling of micro-perforated panel at high sound pressure levels using equivalent fluid approach. *Journal of Sound and Vibration*, 427:134–158, 2018.

[22] Rémi Mignot. *Réalisation en guides d'ondes numériques stables d'un modèle acoustique réaliste pour la simulation en temps-réel d'instruments à vent*. PhD thesis, Télécom ParisTech, 2009.

[23] Gary Paul Scavone. *An acoustic analysis of single-reed woodwind instruments with an emphasis on design and performance issues and digital waveguide modeling techniques*. PhD thesis, 1997.

[24] F. Silva, Ph. Guillemain, J. Kergomard, B. Mallaroni, and A.N. Norris. Approximation formulae for the acoustic radiation impedance of a cylindrical pipe. *Journal of Sound and Vibration*, 322(1-2):255–263, April 2009.

[25] Fabrice Silva, Vincent Debut, Jean Kergomard, Christophe Vergez, Aude Deblevid, and Philippe Guillemain. Simulation of single reed instruments oscillations based on modal decomposition of bore and reed dynamics. *arXiv preprint arXiv:0705.2803*, 2007.

[26] Fabrice Silva, Jean Kergomard, Christophe Vergez, and Joël Gilbert. Interaction of reed and acoustic resonator in clarinetlike systems. *The Journal of the Acoustical Society of America*, 124(5):3284–3295, November 2008. Publisher: Acoustical Society of America.

[27] P-A Taillard, Fabrice Silva, Ph Guillemain, and J Kergomard. Modal analysis of the input impedance of wind instruments. application to the sound synthesis of a clarinet. *Applied Acoustics*, 141:271–280, 2018.

[28] Pierre-André Taillard. *Theoretical and experimental study of the role of the reed in clarinet playing*. PhD thesis, Université du Maine, 2018.

[29] Theodore A Wilson and Gordon S Beavers. Operating modes of the clarinet. *The Journal of the Acoustical Society of America*, 56(2):653–658, 1974.
List of Figures

1 Input impedance of a cylinder $z_{in}$ (computed through Eq. (9) and Eq. (10)) of length $L = 64$ cm.

2 Graphical representation of the $N = 8$ first poles and residues. The values of $v_{RMS}$ are linearly distributed.

3 Error due to the regression of modal coefficients $C_n$ and $s_n$ (Figure 3a), and its consequence on the bifurcation diagram of the $L_2$ norm of $p^{(dim)}$ with respect to $P_m$, for different values of $c_d$, in the absence of a saturation threshold (•).

4 Bifurcation diagram of the $L_2$ norm of $p^{(dim)}$ with respect to $P_m$, for different values of $c_d$, in the presence of a saturation threshold (•).

5 Experimental results from [11] for the saturation threshold (•), compared with three methods:

List of Tables

1 Main constants used for the simulations. 9

2 Values of the extinction threshold $P_{extup} = p_M \cdot \gamma_{extup}$ (kPa), for the simulation of Figure 4a.