A frequency-domain approach of harmonic balance solutions stability
LMA, CNRS, 31 Chemin Joseph Aiguier, 13402 Cedex 20 Marseille, France, e-mail: {vergez,lizee}@lma.cnrs-mrs.fr,
Christophe Vergez, Aude Lize
A method to nd the stability of solutions determined by the harmonic balance method is presented. This method differs from the time-domain Floquet method and takes advantage of intermediate steps of the harmonic balance process. It is exemplied in the case of the self-sustained Van der Pol oscillator.

Introduction

tice to Np harmonics):
Our framework is the study of periodic regimes produced by physical models of self-sustained musical instruments. The harmonic balance method allows to nd both stable and unstable periodic solutions, i.e. to determine through an iterative numerical method their Fourier coefcients as well as the playing frequency. This method is now rather classical in musical acoustics since the work of J. Gilbert et al. ([1]). Moreover a software called Harmbal has been recently developped ([3]) and is freely distributed. However the question of the stability of the solution found by the harmonic balance method is often ignored, though it determines the possibility of observing this solution in real instruments or in time domain numerical methods. When the question of stability is investigated, it is usually done in the time domain, by calculating Floquet multipliers, which can be time-consuming ([2]). Alternatively, in this paper, we propose to use the intermediate results calculated during the harmonic balance process (recalled in section 2). The method relies on a perturbation of this solution like in the classical Floquet approach (section 3.1), but the differential equation fullled by the perturbed solution is written in the frequency domain (see section 3.2). The method is illustrated by studying the well-known self-sustained Van der Pol oscillator, which appears as an extremely simplied model of the clarinet.

n=Np Np

Xn ejnt ,

f (x(t), )

Fn ejnt.
The unknowns to be found by HBM are the Fourier coefcients Xn and the pulsation. Replacing (2) and (3) in (1) leads to : jnXn = Fn XNp ,. , XNp , . Finally, one has to solve G(X) = 0 , (5) (4)
where G is a vector, each component of which being dened n [Np. Np ] by :

Fn XNp ,. , XNp , jnXn.

Harmonic (HBM)

Balance

Method
Practically, the roots of (5) can be found numerically using iteratives techniques. For example in a recent work done by the authors ([3]), the Newton-Raphson scheme is used. This method relies on a linear prediction of the function G, and thus uses its jacobian DXi G, the matrix of rst derivatives of G with respect to the Fourier coefcients Xi. Next section is devoted to the determination of the stability of the periodic solutions found by HBM.
Let a dynamical system be dened by x(t) = f (x(t), ), (1)
Stability of periodic solutions

Time-domain approach

where x is a vector and is the bifurcation parameter. The aim of HBM is to seek periodic solutions x of (1). Since only periodic solutions are considered, terms in (1) can be replaced by their Fourier series (truncated in prac-
Let x (t) be a periodic solution of (1) found by HBM, i.e. a solution of (5). We note xp (t) a perturbation of this solution. The perturbed solution is x(t) x (t) + xp (t) (6)
Forum Acusticum 2005 Budapest

Vergez, Lize

To answer the question of linear stability of x , the nonlinear system (1) is linearized around x. Following Floquets approach ([4] p158), it can be shown that the perturbation may be written as: xp (t) = et (t), (7)
No additional numerical scheme is introduced to answer the question of stability. A similar method has been recently developped in a mechanical context [6]. However, important differences can be noted : In [6], the method is applied to forced oscillations, for which the pulsation is known, on the contrary to self-sustained oscillators. In [6], the method is intended to be applied to systems for which the nonlinear term relies on the stiffness. A direct extension of this method is not straightforward when the nonlinear term includes a resistive part (case of musical instruments, [7]). On the other hand the formalism employed in equation (1) is more general and allows the study of examples with nonlinear resistance (see section 4).

with C, and is a (2/)-periodic function. Rewriting equation (1) for x leads, after linearization, to the fol lowing evolution equation for the perturbation: xp (t) = Dx f |x xp (t) (8)
Replacing xp by its expression (7) then leads to (t) = [ Dx f |x I] (t) (9)
where I is the identity matrix. A practical difculty in (9) is that matrix Dx f |x has periodic coefcients. As explained in ([2]), solving equation (9) requires numerical techniques and artefacts of the numerical schemes introduced may harm the validity of the stability results. Therefore this method is not used. A frequencydomain approach has been developped and is presented below.

Examples

Example : the Van der Pol (VdP) selfsustained oscillator
Frequency-domain approach
The idea is to write the perturbed problem in the Fourier space. Therefore the function and the nonlinear function f are replaced by their Fourier series (truncated to Np harmonics):
The self sustained Van der Pol oscillator is the system : y ry + y = gy 2 y with g > 0, (12)

k ejkt

It is a system particularly interesting in the framework of musical acoustics since it is the archetype of most selfsustained instruments, which may be simplied into oscillators with nonlinear resistance (see e.g. Woodhouse in [7]). 4.1.1 HBM calculus Equation (12) can be reformulated into the form of equation (1) with : x= f= f1 f2 x1 x2 y y , , (13) (14)
Variables of the problem become the Fourier coefcients: det DXk F |X j diagNp kNp (k) I = 0 , k (11) where Xk are Fourier coefcients of the HBM solution x and F is the vector which is constituted by the Fk dened in (3). The unknown in the eigenvalue problem (11) is. If the largest real part of the solutions is strictly positive, then the solution x is unstable, if the largest real part of the solutions is strictly negative, it is stable. The case where the largest real part is zero requires further development of the nonlinear term and is not treated here. Some remarks : A signicant advantage of this method is that matrix DXk F |X j diagN kN (k) is already known k as a by-product of the harmonic balance, since it is the jacobian DXi G as dened in section 2 and evaluated in the last iteration. The use of DXi G to determine the stability of solutions of HBM had already been suggested by J. Gilbert and A. Goncalves ([5]).

x2 gx2 x2 + rx2 x1 1

Let us consider the 1-harmonic approximation of the solution (i.e. Np = 1) and project the nonlinear term onto Fourier basis: F1,1 ejt 1 ejt F1,0 , (15) f1 = F1,1 F2,1 ejt 1 ejt F2,0 , (16) f2 = F2,1

where analytical expression for F1,1 , F1,0 , F1,1 , F2,1 , F2,0 , F2,1 are given in appendix A. In this very simple case, the harmonic balance can be done analytically. Two solutions are exhibited and detailed results are summarized in table 1. Table 1: Harmonic balance results for the 1-harmonic approximation of the VdP oscillator r<0 X1,1 = 0 X1,0 = 0 X1,1 = 0 X2,1 = 0 X2,0 = 0 X2,1 = 0 r0 X1,1 = 0 X1,0 = 0 X1,1 = 0 X2,1 = 0 X2,0 = 0 X2,1 = 0 =1 X1,1 =

X1,0 = 0 X1,1 = r g r g

We note that all the i have a negative real part only for r < 0. Therefore, the non-oscillating solution is stable for r < 0 and unstable for r > 0. The oscillation threshold occurs for r = 0.
solution 1 (non-oscillating)

solution 2 (oscillating)

X1,1 = j X1,0 = 0 X1,1 = j
Figure 1: Non-oscillating solution of the 1-harmonic approximation of the Van der Pol Oscillator: real part of eigenvalues versus damping factor (r in eq. (12))
Oscillating solution In the same way, the evaluation of equation (11) for the oscillating solution gives:
Note that in a typical numerical HBM process as detailed in section 2, the jacobian matrix DXi F (given in appendix B) would be already needed for the calculus of the solution. Therefore it is not an additional cost for the stability analysis. Non-oscillating solution The eigenvalue formulation of the problem (11) for the non-oscillating solution is written as:
Real parts of eigenvalues i are plotted on Figure 2. It can be deduced that the oscillating solution is stable for positive values of the damping factor. This conrms that the stability exchange between the two solutions occurs for r = 0. As expected (Nayfeh, [4], p. 163), one observes that one of the eigenvalues is always zero for the oscillating solution of this automous system (this corresponds to a Floquet multiplier equal to 1). Bifurcation diagram These different results are summarized in the bifurcation diagram presented on Figure 3.
It gives six characteristic multipliers i , whose real parts are plotted on Figure 1 as a function of damping factor r. Let us remark that some of them are degenerated.

0 r+j 0 0

The VdP self-sustained oscillator has been chosen to illustrate the method proposed. The stability results found are therefore well known, and this example should only be considered as a way of validating the method. It is

According to equation (11), the jacobian matrix DXi F dened by equation (25) has to be evaluated for each solution whose stability is under study.

4.1.2 Stability analysis

j 2jr 0 1

0 r + j 0 r

1 = 0. r 0 r j (18)

Conclusion

protable discussions and suggestions.
Expressions of F1,1, F1,0, F1,1, F2,1, F2,0, F2,1
= = = = = = X2,1 X2,0 X2,1 2X1,1 X1,0 X2,0 + 2X1,1 X1,1 X2,1 ) (19) (20) (21) (22) (23) (24)

F1,1 F1,0 F1,1 F2,1

rX2,1 X1,1 g(X1,1 X2,1 + X1,0 X2,1 + 2 rX2,0 X1,0 g(X1,0 X2,0 + 2X1,1 X1,0 X2,1 +

F2,0 F2,1

2X1,1 X1,1 X2,0 + 2X1,0 X1,1 X2,1 )
Figure 2: Oscillating solution of the 1-harmonic approximation of the Van der Pol Oscillator: real part of eigenvalues versus damping factor (r in eq. (12))
rX2,1 X1,1 g(X1,0 X2,1 + X1,1 X2,1 +
2X1,1 X1,1 X2,1 + 2X1,0 X1,1 X2,0 )

Expression of DXi F

DXi F = JA JC JB JD , with (25)

1,1 X2,1

Figure 3: Bifurcation diagram of the 1-harmonic approximation of the Van der Pol oscillator with g = 1, varying the damping factor r. + represents stable equilibrium points whereas o represents unstable equilibrium points. however remarkable that these results are achieved with only one harmonic. Indeed, nothing garantees that the stability of the Np -truncated solution of a system has similar stability properties as the non truncated solution (at least for low Np ). Other information provided by the calculus of the characteristic mupltipliers has not been exploited in this paper. In fact, Im(i ) when max(Re(i )) = 0 allows to identify the type of bifurcation and particularly Hopf or subharmonic bifurcations. In fact, this will be useful when considering more typical models of musical instruments where period-doublings and quasi-periodic regimes are possible. Acknowledgements: the authors want to thank Jol Gilbert, Jean Kergomard and Sergio Bellizzi for very

1 2gX JC

2gX1,0 X2,0 2gX1,1 X2,1 2gX1,0 X2,1 2gX1,1 X2,0
2gX1,0 X2,1 2gX1,1 X2,2gX1,0 X2,0 2gX1,1 X2,1 2gX1,1 X2,1 2gX1,0 X2,1 2gX1,1 X2,0

2gX1,1 X2,1

2gX1,0 X2,1 2gX1,1 X2,0 1 2gX1,1 X2,1

2gX1,1 X2,1 2gX1,0 X2,0

2 gX1,1

2gX1,1 X2,1 2gX1,1 X1,1

JD

2 r gX1,0

2gX1,1 X1,0

2 r gX1,0 2gX1,1 X1,1

2gX1,0 X1,1

2 r 2gX1,0 2gX1,1 X2,1

References
[1] Gilbert J., Kergomard J., and Ngoya E. Calculation of the Steady-State Oscillations of a clarinet using the Harmonic Balance Technique. J. Acoust. Soc. Amer., 86(1):35:41, 1989.
[2] A. Cardona, A. Lerusse, and M. Gradin. Fast fourier nonlinear vibration analysis. Computational Mechanics, 22:128142, 1998. [3] S. Farner, C. Vergez, and J. Kergomard and A. Lize. Contribution to harmonic balance calculations of periodic oscillation for self-sustained musical instruments with focus on single-reed instruments. Submitted to J. Acoust. Soc. Am. Available at http://hal.ccsd.cnrs.fr/ccsd-00004370. [4] A. L. Nayfeh and B. Balachandran. Applied Nonlinear Dynamics: Analytical Computational and Experimental Methods. Wiley series in nonlinear science. Wiley-Interscience, 1995. [5] A. Goncalvs. Stabilit des rgimes doscillations priodiques dinstruments vent. Analyse par quilibrage harmonique. Master Thesis, 2001. [6] G. von Groll and D. J. Ewins. The harmonic balance method with arc-length continuation in rotor/stator contact problems. Journal of Sound and Vibration, 241(2):223233, 2001. [7] J. Woodhouse. Ch5 , Self-Sustained Musical Oscillators. Mechanics of Musical Instruments. Springer Verlag, 185228, 1998.