# Analysis of the self-sustained oscillation of clarinet-like musical instruments using the nonlinear modes approach

Daniel Noreland, Sergio Bellizzi, Robert Bouc, Bruno Cochelin, Jean

Kergomard, Christophe Vergez

► To cite this version:

Daniel Noreland, Sergio Bellizzi, Robert Bouc, Bruno Cochelin, Jean Kergomard, et al.. Analysis of the self-sustained oscillation of clarinet-like musical instruments using the nonlinear modes approach. 2nd International Conference on Nonlinear Normal Modes and Localization in Vibrating Systems, Jun 2006, Samos, Greece. hal-00093638

HAL Id: hal-00093638

https://hal.science/hal-00093638v1

Submitted on 13 Sep 2006

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.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# ANALYSIS OF THE SELF-SUSTAINED OSCILLATION OF CLARINET-LIKE MUSICAL INSTRUMENTS USING THE NONLINEAR MODES APPROACH

D. Noreland, S. Bellizzi, R. Bouc, B. Cochelin, J. Kergomard, C. Vergez

Labratoire de Mécanique et d'Acoustique, CNRS

31 chemin Joseph Aiguier, F-13402 Marseille, France

noreland,bellizzi,bouc,cochelin,kergomard,vergez@lma.cnrs-mrs.fr

## 1 INTRODUCTION

Musical wind instruments are interesting examples of nonlinear vibrating systems. The way in which a self-sustained oscillation is formed is surprisingly complex, even for a simplified model. In short, a wind instrument consists of a resonator and a generator. The resonator is the air column inside the instrument, and is usually characterized by the linear wave equation. The generator consists of some kind of pressure controlled valve, where the relationship between air flow and pressure is starkly non-linear. A simplistic model of sound generation can be based on the assumption that the pressure controlled valve functions as a negative resistance at one of the resonance frequencies of the resonator. Such a model, however, neglects the modal coupling due to the non-linearities that are ultimately the primary reason why different wind instruments sound so different. The aim of this paper is to study how limit cycles of a clarinet-like instrument can be treated in the framework of nonlinear normal modes (NNM). The reason for pursuing this subject is to ultimately be able to derive models of reduced complexity, which are of interest for sound synthesis. Another goal is to identify important control parameters (functions of such entities as blowing pressure, pincing force and position of the player's lips on the reed etc) which can be regulated by a musician in an intuitive way, without a long period of training.

We employ a method devised by Bellizzi and Bouc [11] for computing two-dimensional invariant manifolds of dynamical systems. The choice of this method is dictated by the wish to treat systems with internal resonances. The work presented in this paper should not be seen as a mere application of known methods to a particular example, however. The way the NNMs are treated here is very recent, and it is our hope that part of the work presented will contribute at least to some extent to the knowledge of how to compute the invariant manifolds.

## 2 THE CLARINET MODEL

Wind musical instruments have similar principles of functioning: the player, by blowing inside the instrument destabilizes a valve (a simple reed, a double reed or two lips). The acoustic response of the instrument acts as a feedback loop which influences the valve behavior. The production of a sound corresponds to the auto-oscillation of this dynamical system. Obviously, in spite of these similarities, the functioning of each class of instruments possesses its own specificities. In this section, basic principles of the clarinet functioning are briefly recalled. Simple models have been proposed by [15], [18], [19], [20], or [21] for example of pioneers.

### 2.1 The reed

The reed is often modelled as a mass/spring/damper oscillator. However, because of a resonance frequency ($\approx 3000\$Hz$) large compared to the first harmonics of typical playing frequencies, inertia and damping are often neglected ([18]). This hypothesis leads, considering that reed dynamics is governed by the pressure difference across the reed to:
<figure><img src="image_2.png" /><figcaption>Figure 1: Embouchure of a clarinet (upper part) and its longitudinal cut (lower part).</figcaption></figure>

$$k_s(z - z_0) = (p_{jet} - p_{mouth}) \quad (1)$$

where $z$ ($z_0$) is the reed position (at rest)<fnref n="1" />, $k_s$ is the reed surfacic stiffness, $p_{mouth}$ and $p_{jet}$ are the pressure deviation in the mouth and under the reed tip respectively.

## 2.2 The air flow

As noted by Hirschberg in [22], in the case of clarinet-like instruments, the control of the volume flow by the reed position is due to the existence of a turbulent jet. Indeed, a jet is supposed to form in the embouchure (pressure $p_{jet}$) after the flow separation from the walls, at the end of the (very short) reed channel (see figure 1). Neglecting the velocity of air flow in the mouth compared to jet velocity $v_{jet}$, the Bernoulli theorem applied between the mouth and the reed channel leads to:

$$p_{mouth} = p_{jet} + \frac{1}{2} \rho v_{jet}^2 \quad \text{with } u = S_i v_{jet} \quad (2)$$

where $S_i$ is the cross section at the inlet. Since this area can be expressed as the product between the reed opening $z$ and the reed width $w_r$ (not visible in figure 1 since it is transversal to the plane of the figure), equation (2) can be re-written as follows:

$$u = zw_r \sqrt{\frac{2}{\rho}(p_{mouth} - p_{jet})} \quad (3)$$

where $u$ is the volume flow across the reed, $w_r$ is the width of the reed and $\rho$ is the air density. Combining equations (1) and (3) leads to the wellknown expression of the volume flow as a function of the pressure difference across the reed:

$$u = w_r(z_0 - \frac{1}{k_s}(p_{mouth} - p_{jet}))\sqrt{\frac{2}{\rho}(p_{mouth} - p_{jet})} \quad (4)$$

Since the cross section of the embouchure is large compared to the cross section of the reed channel, it can be supposed that all the kinetic energy of the jet is dissipated through turbulence with no pressure recovery (like in

<footnote n="1">¹Thus the reed is closed when $z = 0$ and opened when $z > 0$</footnote>
the case of a free jet). Therefore, the pressure in the jet is (assuming pressure continuity) the acoustic pressure $p_r$ imposed by the resonator response to the incoming volume flow $u$.

Following [16], we introduce non dimensional pressure $\bar{p} = p_r/(k_s z_0)$ and volume flow $\bar{u} = Z_c u/(k_s z_0)$ and equation (4) writes :

$$ \bar{u} = \zeta(1 + \bar{p} - \gamma)\sqrt{\gamma - \bar{p}} \quad (5) $$

with $\zeta = Z_c l \sqrt{\frac{2z_0}{k_s}}$ and $\gamma = p_{mouth}/(k_s z_0)$. When the reed closes, i.e. $1 + p - \gamma < 0$, the volume is null ($\bar{u} = 0$).

It has been checked in [24] that a cubic expansion of (5) leads to a reasonably good approximation to the resulting periodic solutions, at least below the complete closing of the reed. Therefore, we consider that the volume flow is finally given by :

$$ \bar{u} = u_0 + A\bar{p} + B\bar{p}^2 + C\bar{p}^3 \quad (6) $$

with $u_0 = \zeta(1 - \gamma)\sqrt{\gamma}$, $A = \zeta \frac{3\gamma - 1}{2\sqrt{\gamma}}$, $B = -\zeta \frac{3\gamma + 1}{8\gamma^{2/3}}$, $C = -\zeta \frac{\gamma + 1}{16\gamma^{5/2}}$.

## 2.3 Acoustics in the instrument

We consider a cylindrical bore (length $l$) for the clarinet, and follow [23] p60. The model retained is the wave equation inside the tube, with a source (to take into account the air flow blown inside the instrument) and Neuman and Dirichlet boundary conditions at the input and the output of the tube respectively :

$$ \left\{ \begin{array}{ll} \left[ \partial_{xx}^2 - \left(j \frac{\omega}{c} + \alpha\right)^2 \right] \bar{P}(x) & = -j\omega \frac{\rho}{S} U(\omega) \delta(x), \quad \forall x \in [0; l] \\ \partial_x \bar{P}(x) = 0 & \text{for } x = 0 \\ \bar{P}(x) = 0 & \text{for } x = l \end{array} \right. \quad (7) $$

The dimensionless pressure field $\bar{p}(x, t)$ is decomposed onto the modes of the air column inside the bore :

$$ \bar{p}(x, t) = \sum_{n=1}^{\infty} f_n(x) p_n(t) \quad (8) $$

where the family {$f_n$}$_{n \in \mathbb{N}}$ is a basis of orthogonal eigenmodes of air column. In the case of a close/opened cylindrical bore with dispersion neglected $f_n(x) = \cos k_n(x)$, where $k_n = \frac{n\pi}{2l}$ and $n$ is an odd positive integer and $l$ the length of the bore.

Modal coordinates $p_n(t)$ are calculated through the projection of equation (7) written in the time domain (with $\bar{p}$ replaced by its expansion (8) truncated to $N$ modes) on each mode {$f_n$}.

$$ \ddot{p}_n(t) + 2\alpha c \dot{p}_n(t) + (\omega_n^2 - \alpha^2 c^2) p_n(t) = \frac{2c}{l} \dot{u}(t) \quad (9) $$

where $\omega_n = k_n c$. As explained p70 in [23], the damping of each mode is low ($\alpha c \ll \omega_n$) so that the third term of the left member can be reduced. Moreover, following [23] p72, we consider that the value of $\alpha$ associated to each mode may be different and is noted $\alpha_n \simeq \frac{Y_n}{l}$ where $Y_n$ is the value of the admittance at frequency $\omega_n/2\pi$. Finally, equation (9) is rewritten :

$$ \ddot{p}_n(t) + 2Y_n \frac{c}{l} \dot{p}_n(t) + \omega_n^2 p_n(t) = \frac{2c}{l} \dot{u}(t) \quad (10) $$

## 2.4 Complete model

Considering equations (10) and (6) leads to the dimensionless model made of $N$ second order ODE :

$$ \ddot{p}_n(t) + 2Y_n \frac{c}{l} \dot{p}_n(t) + \omega_n^2 p_n(t) = \frac{2c}{l} \left( A + 2B \sum_{i=1}^{N} p_i(t) + 3C \left( \sum_{i=1}^{N} p_i(t) \right)^2 \right) \sum_{i=1}^{N} \dot{p}_i(t) \quad (11) $$
with $A = \zeta \frac{3\gamma - 1}{2\sqrt{\gamma}}$, $B = -\zeta \frac{3\gamma + 1}{8\gamma^{2/3}}$, $C = -\zeta \frac{\gamma + 1}{16\gamma^{5/2}}$

The intention now is to apply the concept of nonlinear modes to this equation.
## 3 NONLINEAR MODES

We briefly recall here how to characterize the Nonlinear Modes (NMs) in the framework of the invariant manifold theory using an amplitude-phase transformation. A complete presentation can be found in [8][9][10].

We consider a system of the form

$$[M] \ddot{Q}(t) + F(\dot{Q}(t), Q(t)) = 0 \quad (12)$$

where $Q(t)$ is an $n$-vector, $[M]$ is a non-singular symmetric square $n \times n$-matrix and $F$ is a (sufficiently regular) vector function with dimension $n$ such that $F(0, 0) = 0$.

### 3.1 Definition of a nonlinear mode

We focus on motions (solutions of Eq. (12)) where the displacements and velocities, $(Q(t), \dot{Q}(t))$ are related to a single pair of amplitude and phase variables, $(v(t), \phi(t))$, according to

$$\left\{ \begin{aligned} \dot{Q}(t) &= v(t)\mathbf{X}(v(t), \phi(t)) \\ \dot{\phi}(t) &= \Omega(v(t), \phi(t)) \\ \dot{v}(t) &= v(t)\xi(v(t), \phi(t)) \end{aligned} \right. \quad \text{with } \phi(0) = \varphi \text{ and } v(0) = a \quad (13)$$

where $\mathbf{X}$ and $\mathbf{Y}$ are $n$-vector functions, which are $2\pi$-periodic with respect to the phase variable. The amplitude and phase variables are governed by the two first order differential equations

$$\left\{ \begin{aligned} \dot{\phi}(t) &= \Omega(v(t), \phi(t)) \\ \dot{v}(t) &= v(t)\xi(v(t), \phi(t)) \end{aligned} \right. \quad \text{with } \phi(0) = \varphi \text{ and } v(0) = a \quad (14)$$

where $\Omega$ (the frequency function or frequency modulation function) and $\xi$ (the damping function or amplitude modulation function) are two scalar functions, which are even and $\pi$-periodic with respect to the phase variable and $\varphi \in (0, 2\pi]$ and $a > 0$ are two given constants which set the initial conditions of the motion.

If such motions (13)(14) exist, they define a nonlinear mode for Eq. (12) which is characterized by the four functions $\mathbf{X}$, $\mathbf{Y}$, $\Omega$ and $\xi$. The modal motion is confined to lie on a two-dimensional invariant manifold ([4][13]) in the phase space, defined by the parametric equations

$$\left\{ \begin{aligned} \dot{\mathbf{Q}} &= a\mathbf{X}(a, \varphi) \\ \dot{\mathbf{\phi}} &= a\mathbf{Y}(a, \varphi) \end{aligned} \right. \quad (15)$$

The dynamics of the motion on the invariant manifold is given by Eqs. (14). If $\xi \neq 0$ (i.e. $v(t) = a, \forall t$), the modal motions defined by Eqs. (13) and (14) will be periodic with a period given by

$$T(a) = \int_{0}^{2\pi} \frac{1}{\Omega(a, \phi)} d\phi. \quad (16)$$

Periodic solutions may also exist if $\xi/\Omega \neq 0$ does not keep a constant sign. Indeed, from Eq.(13), it follows that

$$\frac{dv}{d\phi} = v\tau(v, \phi) \quad (17)$$

where $\tau(v, \phi) = \xi(v, \phi)/\Omega(v, \phi)$ can be viewed as a "generalized damping rate function". Since the right hand side of Eq. (17) is $\pi$-periodic with respect to $\phi$, periodic solutions $(v^*(\phi) = v^*(\phi + \pi))$ may exist with some $\xi$ and $\Omega$ (one necessary condition being that $\tau(v, \phi)$ does not keep a constant sign). It follows that

$$\mathbf{Q}(t) = v^*(\phi(t))\mathbf{X}(v^*(\phi(t)), \phi(t))$$

will be a $T$-periodic function with period

$$T = \int_{0}^{2\pi} \frac{d\phi}{\Omega(v^*(\phi), \phi)}, \quad (18)$$

and the stability analysis of this periodic function (i.e. a "limit cycle" when Eqs. (13)(14) are a solution in the case of a nonlinear oscillator) will be easily deduced from the variational equation associated with Eq. (12).
## 3.2 Characterization of a nonlinear mode

Substituting Eq. (13) into Eq. (12) and using Eq. (14) and the chain rule yield a set of first order nonlinear Partial Differential Equations (PDEs) in the two variables $(v, \phi)$,

$$ (\mathbf{X} + v\mathbf{X}_v) \xi + \mathbf{X}_\phi \Omega = \mathbf{Y} \quad (19) $$

$$ [\mathbf{M}] (\mathbf{Y} + v\mathbf{Y}_v) \xi + [\mathbf{M}] \mathbf{Y}_\phi \Omega + \frac{1}{v} \mathbf{F}(v\mathbf{Y}, v\mathbf{X}) = \mathbf{0}. \quad (20) $$

where $(.)_\phi$ denotes the partial differentiation with respect to $\phi$. The PDEs (19)(20) are independent of time and contain the unknown vector functions $\mathbf{X}$ and $\mathbf{Y}$ its partial derivatives with respect to $v$ and $\phi$ and the two unknown scalar functions $\Omega$ and $\xi$.

In order to characterize the four unknown functions $\mathbf{X}$, $\mathbf{Y}$, $\Omega$ and $\xi$, it is necessary to add two constraint scalar equations to (19) and (20) (often called normalization conditions). Due to the $2\pi$-periodicity with respect the variable $\phi$, the functions $\mathbf{X}$ can be expressed as

$$ \mathbf{X} = \mathbf{X}^{oc} + \mathbf{X}^{ec} + \mathbf{X}^{os} + \mathbf{X}^{es} \quad (21) $$

where $(.)^{oc}$, $(.)^{ec}$, $(.)^{os}$, $(.)^{es}$, respectively) denotes the odd cosine (even cosine, odd sine, even sine, respectively) terms in the corresponding Fourier series. We will adopt in this study, without any loss of generality, the following constraint scalar equations

$$ \mathbf{X}^{ocT} [M] \mathbf{X}^{oc} + \mathbf{X}_\phi^{osT} [M] \mathbf{X}_\phi^{os} = \cos^2 \phi, \quad (22) $$

$$ \mathbf{X}^{ocT} [M] \mathbf{X}_\phi^{os} = 0. \quad (23) $$

These constraint equations reduce to usual normalization conditions for some important special cases, for example the linear case in the next section.

A nonlinear mode of the system (12) is obtained by solving Eqs. (19)(20)(22)(23) for the four functions $\mathbf{X}$, $\mathbf{Y}$, $\Omega$ and $\xi$, recalling that

$$ \mathbf{X}(v, \phi) = \mathbf{X}(v, \phi + 2\pi), \quad \Omega(v, \phi) = \Omega(v, -\phi) = \Omega(v, \phi + \pi) \quad \text{and} \quad \xi(v, \phi) = \xi(v, -\phi) = \xi(v, \phi + \pi). \quad (24) $$

It is worth noting that depending on the properties of function $\mathbf{F}$, some functions among $\mathbf{X}^{os}$, $\mathbf{X}^{es}$, $\mathbf{X}^{oc}$, $\mathbf{X}^{ec}$ can be discarded. For example,

$$ \begin{array}{l@{\quad}l@{\quad}l} \text{if } \mathbf{F}(\mathbf{Q}, \dot{\mathbf{Q}}) = \mathbf{F}(\mathbf{Q}) & \text{then} & \mathbf{X}^{os} \equiv \mathbf{0}, \quad \mathbf{X}^{es} \equiv \mathbf{0}, \\ \text{if } \mathbf{F}(\mathbf{Q}, \dot{\mathbf{Q}}) = -\mathbf{F}(-\mathbf{Q}, -\dot{\mathbf{Q}}) & \text{then} & \mathbf{X}^{ec} \equiv \mathbf{0}, \quad \mathbf{X}^{es} \equiv \mathbf{0}. \end{array} $$

## 3.3 Linear case

If $\mathbf{F}(\mathbf{Q}, \dot{\mathbf{Q}}) = [\mathbf{K}] \mathbf{Q} + [\mathbf{C}] \dot{\mathbf{Q}}$ where $[\mathbf{K}]$ and $[\mathbf{C}]$ are square matrices, Eqs. (19)(20)(22)(23) can be easily solved given

$$ \left\{ \begin{array}{l} \mathbf{X}(v, \phi) = \psi^c \cos \phi - \psi^s \sin \phi \\ \Omega(v, \phi) = \omega, \text{ and } \xi(v, \phi) = \eta \end{array} \right. \quad (25) $$

where $\omega$, $\eta$, $\psi^c$ and $\psi^s$ satisfy the eigenvalue problem

$$ \left( \begin{array}{cc} [\mathbf{C}] & [\mathbf{M}] \\ [\mathbf{M}] & [\mathbf{0}] \end{array} \right) \mathbf{\Psi} \lambda + \left( \begin{array}{cc} [\mathbf{K}] & [\mathbf{0}] \\ \mathbf{0} & -[\mathbf{M}] \end{array} \right) \mathbf{\Psi} = \mathbf{0} $$

with $\mathbf{\Psi} = (\psi^T, \lambda\psi^T)^T$, $\lambda = \eta + i\omega$ (assuming $\omega \neq 0$) and $\psi = \psi^c + i\psi^s$ with the normalization condition

$$ \begin{array}{l} \psi^{cT} [M] \psi^c + \psi^{sT} [M] \psi^s = 1, \\ \psi^{cT} [M] \psi^s = 0. \end{array} $$

Hence the formulation (13)(14) can be viewed as an extension of the linear modal analysis.
# 4 Numerical solution of the equations describing the manifold

Eq. (19), (20), (22), and (23) constitute a partial differential algebraic equation (PDAE), where $v$ acts as a time-like variable, and the problem has periodic boundary conditions in the $\phi$-direction. One way of solving this PDAE is to discretise it in the $v$- and $\phi$-directions, and to employ finite difference approximations of the derivatives. This approach differs from those in [7] or [10], where use is made of Galerkin methods based on trigonometric terms, and in the latter case polynomial terms in the $v$-direction. Although elegant, the Galerkin treatment becomes prohibitively complex as the number of expansion terms increases. A low order, implicit scheme for step-wise advancement in the $v$-direction was also deemed to be better adapted for capturing irregularities in the solution, which do not readily lend themselves to an accurate description with a polynomial basis.

Since the problem is periodic in the $\phi$-direction, and can be expected to have a smooth solution in this direction, it is natural to use pseudo-spectral approximations[14] of the $\partial/\partial\phi$-terms. This pseudo-spectral approximation is closely related to a Galerkin method with the usual trigonometric basis, but can also be thought of as a usual finite difference approximation where the number of points, and hence the order of accuracy, is a function of the step-size $h_\phi$ in the direction of differentiation. For sufficiently smooth functions, the resulting approximation has spectral convergence rate, meaning that the error decreases faster than any polynomial in $h_\phi$.

The equation is left stated on its implicit form, since it may have physically interesting solutions for which some $\xi(v, \phi) = 0$. By approximating the derivatives in eq. (19), (20), (22), and (23), with the pseudo-spectral scheme in the $\phi$-direction and a backward difference in the $v$-direction, an implicit Euler approximation of the PDAE at hand is obtained. The solution on each new $v$-level is then found numerically using the Newton method.

For the solutions considered in this work, the non-linear modes are the continuations of corresponding linear modes. Consistent initial conditions of the PDAE are therefore easily obtained from the linear system, or simply by setting $v = 0$ in eq. (19)-(23), which then collapses to an algebraic system of equations.

## 4.1 Computation of time dependent solutions from the manifolds

Knowing the manifolds $X(v, \phi)$, $Y(v, \phi)$ and the damping and frequency functions $\xi(v, \phi)$ and $\Omega(v, \phi)$, the time evolutions $v(t)$ and $\phi(t)$ can be computed numerically solving eq. (14). Since numerical approximations of $\xi$ and $\Omega$ are known only for certain values $v_i$ and $\phi_i$, a two-dimensional interpolation procedure is used, employing trigonometric interpolation in the $\phi$-direction and quadratic interpolation in the $v$-direction. For solutions representing a limit cycle, it is necessary to solve eq. (19)-(23) on an interval $[0, v_{max}]$, where $v_{max}$ is slightly above the amplitude of the limit cycle. An estimate of this limit can be obtained by keeping track of the mean damping function

$$ \langle \tau(v) \rangle = \frac{1}{2\pi} \int_{0}^{2\pi} \frac{\xi(v, \phi)}{\Omega(v, \phi)} d\phi, $$

which is zero near the oscillation threshold. The physical variables in phase space are finally given by eq. (13), where once again the two-dimensional interpolation procedure is employed.

## 4.2 The numerical scheme applied on the clarinet model

We consider the case of a clarinet model as described in section 2, with three modes (N=3). The method presented in section 3 is applied to find the NLMs of the system. This corresponds to $Q = [p_1 p_2 p_3]^T$. The model parameters are chosen so as to correspond to a mf playing condition. In a first step, the invariant manifolds describing the components of $p$ and $\dot{p}$, as well as $\xi$ and $\Omega$, are computed. The shapes of the manifolds are shown in Fig. (2) for a case with 31 discretization points in the $\phi$-direction (effectively resolving the first 15 Fourier terms), and 50 discretization points on the interval $v \in [0, 0.5]$. As expected, the shape of $p_1$ starts out from a purely trigonometric shape at $v = 0$. The shape then changes only slightly as the amplitude grows. Components $p_2$ and $p_3$ are zero at $v = 0$, but then change dramatically as the amplitude grows.

In a second step, the time evolutions $v(t)$ and $\phi(t)$ are calculated, and from the results, also the time evolutions of the components of $p(t)$ and $\dot{p}(t)$. These results are shown in Fig. (4)-(5). The simulation of $v$ and $\phi$ has been started from a rather arbitrary, non-zero initial condition $v(0) = 0.1$, $\phi(0) = 0$. As can be seen, the amplitude of the
<figure><img src="image_3.png" /><figcaption>Figure 2: The two-dimensional manifolds for the three components of X and Y. Shown at the bottom are also the damping and amplitude modulation functions.</figcaption></figure>

first component grows quickly initially, but then stabilize as the limit cycle is approached. The same phenomenon is visible for the higher modes, but they show relatively a much stronger growth. This is a typical feature for wind instruments, where the small amplitude oscillations are nearly sinusoidal, whereas nonlinear effects add increasingly to the timbre as the amplitude grows. In terms of the manifolds, this is seen both as increasing amplitude of the higher components ($n > 1$) and a minute variation from the harmonic shape.

In order to check the feasibility of the solutions, a reference solution was computed by solving Eq. (11) in Matlab using `ode45` with a small tolerance. The difference between the solution obtained from the manifolds, and the reference solution is small, and can be seen in Fig. (6)–(7). As can be expected, the error grows slowly in time due to the numerical phase error in the solution.

## 4.3 Future work

The work presented in this paper is still in progress. Future work entails the investigation of well-posedness of the continuous equations (19)–(23). Some problems with divergence of the numerical solution has been observed for certain values of $v$, but it remains to investigate whether this is a numerical problem, or a sign of ill-posedness of the underlying equations. There are indications that the value $v_i$, where the instability occurs, converges to a certain value as $h_v$ decreases, indicating a feature of the continuous equation, rather than a numerical problem. Another thing that remains to be investigated is how bifurcations of NNMs affect the solution.

As far as the numerical integration procedure is concerned, better use of the flexibility of the one-step method used can be imagined, involving i.e automatic step-size control based on a local error estimates.
<figure><img src="image_4.png" /><figcaption>Figure 3: Time evolution $v(t)$ and $\phi(t)$. Not clearly visible in the figures is a fine, periodic ripple in $v$ and $\phi$. The ripple has the same period as the limit cycle. The bottom curve shows the mean damping function $\langle \tau(v) \rangle$, whose zero indicates the amplitude of the limit cycle.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: The solutions of the three components $p_1, p_2$ and $p_3$.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: The solutions of the three components $\dot{p}_1, \dot{p}_2$ and $\dot{p}_3$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: The absolute errors of the three components $p_1, p_2$ and $p_3$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 7: The absolute errors of the three components $\dot{p}_1$, $\dot{p}_2$ and $\dot{p}_3$.</figcaption></figure>

References

[1] R.M. Rosenberg, On non linear vibrations of systems with many degree of freedom. *Advances in Applied Mechanics*, **242**(9), 155-242, 1966.

[2] W. Szemplinska-Stupnicka, *The Behavior of Nonlinear Vibrating Systems* (Volumes I and II). Kluwer Academic Publishers, 1990.

[3] L. Jezequel, C.H. Lamarque, Analysis of nonlinear dynamical systems by normal form theory. *Journal of Sound and Vibration*, **149**(3), 429-459, 1991.

[4] S.W. Shaw and C. Pierre, Non-linear normal modes and invariant manifolds. *Journal of Sound and Vibration*, **150**(1), 170-173, 1991.

[5] A.F. Vakakis Nonlinear normal modes (NNMs) and their applications in vibration theory. *Mechanical Systems and Signal Processing*, **11**(1), 3-22, 1997.

[6] S.W. Shaw and C. Pierre, Normal modes for nonlinear vibratory systems. *Journal of Sound and Vibration*, **164**(1), 85-124, 1993.

[7] E. Pescheck, S.W. Shaw and C. Pierre, A new Galerkin-based approach for accurate nonlinear normal modes through invariant manifolds. *Journal of Sound and Vibration*, **249**(5), 971-993, 2002.

[8] S. Bellizzi and R. Bouc, A new formulation for the existence and calculation of non-linear normal modes. *Journal of Sound and Vibration*, **287**, 545-569, 2005.

[9] S. Bellizzi and R. Bouc, A New Formulation for the Existence and Calculation of Nonlinear Modes. S. Bellizzi, B. Cochelin, C.H. Lamarque eds. *EUROMECH 457 Nonlinear Modes of Vibrating Systems, Press of Ecole Nationale des Travaux Publics de l'Etat* (ISBN 2-86 834-119-5), 25-28, 2004.

[10] S. Bellizzi and R. Bouc, An amplitude phase formulation for nonlinear modes and limit cycles through invariant manifolds. *Journal of Sound and Vibration*, (submitted) 2006.
[11] S. Bellizzi and R. Bouc, Nonlinear modes: amplitude-phase formulation and limit cycles analysis. *2nd Int. Conf. on Nonlinear Normal Modes and Localization in Vibrating Systems, Samos, June 19-23, 2006*.

[12] J. Hale, *Ordinary Differential Equations*. Wiley-Interscience, New York, 1969.

[13] R. Arquier, F. Pérignon and B. Cochelin, Numerical Continuation of Non Linear Modes of Elastic Structures. S. Bellizzi, B. Cochelin, C.H. Lamarque eds. *EUROMECH 457 Nonlinear Modes of Vibrating Systems, Press of Ecole Nationale des Travaux Publics de l'Etat* (ISBN 2-86 834-119-5), 17-20, 2004.

[14] B. Gustafsson, H. O. Kreiss and J. Oliger, *Time dependent problems and difference methods*. John Wiley & Sons, New York, 1995.

[15] J. Backus. Small vibration theory of the clarinet. *J. Acoust. Soc. Amer.*, **35**, 305, 1963. and Erratum (61) [1977], 1381.

[16] J. Kergomard. *Mechanics of Musical Instruments (Chap. 6)*, chapter Elementary considerations on reed-instrument oscillations. Springer Verlag, 1995.

[17] A. Hirschberg. *Mechanics of Musical Instruments (Chap. 7)*, chapter Aero-acoustics of wind instruments. Springer Verlag, 1995.

[18] C. J. Nederveen. *Acoustical aspects of woodwind instruments*. Fritz Knuf pub., Amsterdam, 1969.

[19] W. E. Worman. *Self-sustained nonlinear oscillations of medium amplitude in clarinet-like systems*. PhD thesis, Case Western Reserve University, 1971. Ann Arbor University Microfilms (ref. 71-22869).

[20] T. A. Wilson and G. S. Beavers. Operating modes of the clarinet. *J. Acoust. Soc. Amer.*, **56**, 653, 1974.

[21] R. T. Schumacher. Ab initio calculations of the oscillations of a clarinet. *Acustica*, 48:71-85, 1981.

[22] A. Hirschberg, J. Gilbert, A. P. J. Wijnands, and A. M. C. Valkering. Musical aero-acoustics of the clarinet. *Journal de Physique IV*, 4:C5-559:C5-568, Mai 1994. Colloque C5 supplément au Journal de Physique III.

[23] V. Debut, Deux études d’un instrument de musique de type clarinette : analyse des fréquences propres du résonateur et calcul des auto-oscillations par décomposition modale, phd (dir. J. Kergomard and F. Laloë). Aix-Marseille 2 university.

[24] C. Fritz and S. Farner and J. Kergomard, Some aspects of the harmonic balance method applied to the clarinet, *Applied acoustics*, 65:1155-1180, 2004