n-dimensional Harmonic Balance Method extended to
non-explicit nonlinearities

Mathias Legrand, Sébastien Roques, Christophe Pierre, Patrice Cartraud,
Bernard Peseux

► To cite this version:

Mathias Legrand, Sébastien Roques, Christophe Pierre, Patrice Cartraud, Bernard Peseux. n-dimensional Harmonic Balance Method extended to non-explicit nonlinearities. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2006, 15 (1-2-3), pp.269-280. 10.3166/remn.15.269-280. [https://hal.science/hal-00354486v2](https://hal.science/hal-00354486v2)

HAL Id: hal-00354486

https://hal.science/hal-00354486v2

Submitted on 28 Apr 2009

**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
# *n*-dimensional Harmonic Balance Method extended to
non-explicit nonlinearities

**Mathias Legrand**<fnref n="1" />

GÉM - Institut de Recherches en Génie Civil et Mécanique UMR CNRS 6183, équipe *Structure et Simulation*, École Centrale de Nantes, 1 rue de la Noë 44321 Nantes Cedex 3, France

**Sébastien Roques**

GÉM - Institut de Recherches en Génie Civil et Mécanique UMR CNRS 6183, équipe *Structure et Simulation*, École Centrale de Nantes, 1 rue de la Noë 44321 Nantes Cedex 3, France

**Christophe Pierre**

Structural Dynamics and Vibration Laboratory, Department of Mechanical Engineering,
McGill University, 817 Sherbrooke St West, Montreal, Quebec H3A 2K6, Canada

**Patrice Cartraud**

GÉM - Institut de Recherches en Génie Civil et Mécanique UMR CNRS 6183, équipe *Structure et Simulation*, École Centrale de Nantes, 1 rue de la Noë 44321 Nantes Cedex 3, France

**Bernard Peseux**

GÉM - Institut de Recherches en Génie Civil et Mécanique UMR CNRS 6183, équipe *Structure et Simulation*, École Centrale de Nantes, 1 rue de la Noë 44321 Nantes Cedex 3, France

**Abstract**

The harmonic balance method is widely used for the analysis of strongly nonlinear prob-
lems under periodic excitation. The concept of *hypertime* allows for the generalization
of the usual formulation to multi-tone excitations. In this article, the method is applied
to a system involving a nonlinearity which cannot be explicitly expressed in the multi-
frequency domain in terms of harmonic coefficients. The direct and inverse Discrete Fast
Fourier Transforms are then necessary to alternate between time and frequency domains
in order to take into account this nonlinearity. The results show the efficiency and the pre-
cision of the method.

**Keywords:** harmonic balance method, hypertime domain, unilateral contact

Résumé

La méthode de l'équilibrage harmonique permet l'analyse de problèmes fortement non
linéaires sous excitation périodique. En utilisant le concept d'hyper-temps, il est possible
de généraliser la formulation usuelle à des excitations apériodiques. Dans cet article, cette
méthode est appliquée à un système impliquant une non-linéarité qui ne s'exprime pas de
façon explicite dans l'espace multi-fréquentiel des coefficients harmoniques. Les transfor-
mées de Fourier rapides directe et inverse sont alors nécessaires pour alterner domaines
temporel et fréquentiel afin de prendre en compte cette non-linéarité. Les résultats mon-
trent l'efficacité et la précision de la méthode.

**Mots clés :** méthode de l'équilibrage harmonique, espace hyper-temps, contact uni-
latéral

<footnote n="1">¹Corresponding author: mathias.legrand@ec-nantes.fr</footnote>
# 1 Introduction

It is well known the the Harmonic Balance Method (HBM) is one of the most commonly used approaches for analyzing strong nonlinear systems. It offers an alternative to time-domain methods for the analysis of equations where a periodic steady-state is sought and has been used for years in electrical engineering for studying nonlinear circuits [1]. The usual one-dimensional HBM is limited to predicting periodic vibrations only whereas aperiodic (quasi-periodic) solutions made up with incommensurate tones are frequent. For the first time in [2], the HBM applied to mechanical systems was extended to the multi-tone free response of a nonlinear beam with cubic terms considering the concept of hypertime. Unfortunately, the main drawback of this approach, which requires a complicated algebraic formulation as a preliminary, lies in its inability to consider nonlinearities, such as unilateral contact, that can not be directly evaluated in the frequency domain. A major breakthrough was proposed by [3] as the Alternating Frequency/Time domain method (AFT). It was realized that strong nonlinearities can be accurately analyzed in the time domain and transformed back in the frequency domain so as to form a set of nonlinear equations which can be handled by a nonlinear solver. Based on a pioneering work [4], an extension of the AFT to aperiodic systems was proposed in [5] for the study of the internal resonant vibration of a nonlinear Jeffcott rotor with contact terms using multi-frequency Fourier Transforms. In a recent work, Pušenjak [6] has extended the method to handle general multi-degree of freedom externally excited and autonomous dynamical systems.

This article finds its motivation in a study concerning the rotor-stator "modal" interaction explained in [7]. According to this work, for a tiny range of rotational velocities and under certain circumstances, a travelling wave speed coincidence [8] will occur between the nodal diameter modes of the bladed-disk and the casing through direct contact, leading to an aperiodic limit cycle. Time integration including inherent difficulties like time-step size and initial conditions, it has been found interesting to use the presented n-dimensional HBM for a clearer insight and an intrinsic description of this phenomenon.

# 2 n-dimensional harmonic balance method: theoretical statements

Nonlinear mechanical problems generally involve equations of motion which can be written in the following generic form:

$$M\ddot{X}+D\dot{X}+KX+F_{nl}(X,\dot{X})=F_{ext} \quad (1)$$

where X, $\dot{X}$ and $\ddot{X}$ respectively stand for the position, velocity and the acceleration vectors expressed in generalized coordinates. The external excitation is given in $F_{ext}$ and the nonlinear terms are gathered in $F_{nl}(X, \dot{X})$. The mass, damping and stiffness matrices of the mechanical system are respectively named M, D and K. Following [6], the extension of the HBM to multi-tone systems requires the construction of an hypertime domain. First, $M_s$ pseudo-time variables $\tau_m = \omega_m t$ for $m = 1, 2, \dots, M_s$ within the hypertime domain $0 \le \tau_m \le 2\pi$ are introduced, where $\omega_m$ represent the incommensurable frequencies of the steady-state response sought. Assuming that the aperiodic solution for the generalized coordinate $X_j$ (jth coordinate of X for $j = 1, 2, \dots, N$) exists where N represents the total number of DOF of the system, the solution is sought in the form of truncated multiple Fourier series:

$$X = \sum_{j_1=-N_h}^{N_h} \sum_{j_2=-N_h}^{N_h} \dots \sum_{j_m=-N_h}^{N_h} \left[ a_{j_1, j_2, \dots, j_{M_s}} \cos \left( \sum_{m=1}^{M_s} j_m \tau_m \right) + b_{j_1, j_2, \dots, j_{M_s}} \sin \left( \sum_{m=1}^{M_s} j_m \tau_m \right) \right] \quad (2)$$

where $N_h$ denotes the highest order of the harmonic terms in the truncated Fourier series. In equation (2), the various combination tones are expressed by means of the arguments $\sum_{m=1}^{M_s} j_m \tau_m$,
where $j_m$ are positive or negative integers. Vectors $\mathbf{a}_{j_1, j_2, \dots, j_{M_s}}$ and $\mathbf{b}_{j_1, j_2, \dots, j_{M_s}}$ include all the Fourier coefficients of the frequency-discretized displacement along with the different tone combinations. Because of the parity of the cos and sin functions, only the terms at positive combination frequencies meeting the following constraint need to be kept:

$$ \sum_{m=1}^{M_s} j_m \tau_m \ge 0 \qquad (3) $$

In spite of the above restriction on the positive combination frequencies, the total number of harmonic terms in the solution will become very high at greater values of $N_h$. To keep the smallest total number of harmonic terms, additional constraints on the combination frequencies are assumed if possible. Once the selection of combination frequencies and harmonics has been performed, the solution can be constructed by following equation (2) which makes possible to express the aperiodic displacement $X_j$ conveniently in a matrix form as:

$$ X_j = \mathbf{T} \cdot \mathbf{u}_j, \quad j = 1, 2, \dots, N \qquad (4) $$

where $\mathbf{T} = [\mathbf{T}^{\cos}, \mathbf{T}^{\sin}]$, $\mathbf{u}_j = \{\mathbf{u}_{\cos}^j, \mathbf{u}_{\sin}^j\}^T$ and superscript $^T$ denotes matrix transposition. At this level, it is worth mentioning that $\mathbf{u}_j$ is a non-redundant rearrangement of the harmonic coefficients $a_{j_1, j_2, \dots, j_{M_s}}$ and $b_{j_1, j_2, \dots, j_{M_s}}$ defined in equation (2) and consistent with the row vector $\mathbf{T}$, combination of two subvectors $\mathbf{T}^{\cos}$ and $\mathbf{T}^{\sin}$ of the form:

$$ \mathbf{T}^{\cos} = \left[ \cos\left(\sum_{m=1}^{M_s} j_m^1 \tau_m\right), \cos\left(\sum_{m=1}^{M_s} j_m^2 \tau_m\right), \dots, \cos\left(\sum_{m=1}^{M_s} j_m^{N_{\cos}} \tau_m\right) \right] \qquad (5) $$

where $j_m^i$ ($i = 1, 2, \dots, N_{\cos}$ and $m = 1, 2, \dots, M_s$) are integers satisfying constraint (3). Subvector $\mathbf{T}^{\sin}$ is analogously combined of $N_{\sin}$ elements of sine terms. The total number of harmonic terms $N_h$ is equal to $N_{\cos} + N_{\sin}$. Considering that the constant component can be interpreted by the corresponding cosine term at the zero frequency in matrix $\mathbf{T}$ and that no sine term at the zero frequency in matrix $\mathbf{T}$ is needed yields: $N_{\cos} = N_{\sin} + 1$. By introducing the Kronecker product $\otimes$ such as:

$$ \mathbf{Y} = \mathbf{T} \otimes \mathbf{I}_N \quad \text{and} \quad \mathbf{u} = [\mathbf{u}_1^T, \mathbf{u}_2^T, \dots, \mathbf{u}_N^T]^T \qquad (6) $$

vectors $\mathbf{X}$, $\dot{\mathbf{X}}$ and $\ddot{\mathbf{X}}$ are expressed as follows:

$$ \mathbf{X} = \mathbf{Y} \cdot \mathbf{u}, \quad \dot{\mathbf{X}} = \left( \sum_{n=1}^{M_s} \omega_n \frac{\partial \mathbf{Y}}{\partial \tau_n} \right) \cdot \mathbf{u} \quad \text{and} \quad \ddot{\mathbf{X}} = \left( \sum_{m=1}^{M_s} \sum_{n=1}^{M_s} \omega_m \omega_n \frac{\partial^2 \mathbf{Y}}{\partial \tau_m \partial \tau_n} \right) \cdot \mathbf{u} \qquad (7) $$

and plugged into equation (1) before a Galerkin procedure is applied over the entire hyper-time domain $[0, 2\pi]^{M_s}$, resulting in a vector of nonlinear equations in $\mathbf{u}$ that are to be solved simultaneously:

$$ \mathbf{R} = \mathbf{H} \cdot \mathbf{u} + \mathbf{F}_{nl}^{u}(\mathbf{u}) - \mathbf{F}_{ext}^{u} \qquad (8) $$

Notations in equation (8) are given below:

$$ \mathbf{H} = \left\langle \mathbf{Y}, \left[ \sum_{m=1}^{M_s} \omega_m \left( \sum_{n=1}^{M_s} \omega_n \mathbf{M} \frac{\partial^2 \mathbf{Y}}{\partial \tau_m \partial \tau_n} + \mathbf{D} \frac{\partial \mathbf{Y}}{\partial \tau_m} \right) + \mathbf{K} \mathbf{Y} \right] \right\rangle \qquad (9) $$

$$ \mathbf{F}_{nl}^{u}(\mathbf{u}) = \langle \mathbf{Y}, \mathbf{F}_{nl}(\mathbf{X}, \dot{\mathbf{X}}) \rangle $$

$$ \mathbf{F}_{ext}^{u} = \langle \mathbf{Y}, \mathbf{F}_{ext} \rangle $$

with the inner product $\langle \cdot, \cdot \rangle$ defined as follows:

$$ \langle \mathbf{Y}_1, \mathbf{Y}_2 \rangle = \frac{2}{(2\pi)^{M_s}} \int_0^{2\pi} \int_0^{2\pi} \dots \int_0^{2\pi} \mathbf{Y}_1^T \mathbf{Y}_2 \, d\tau_1 d\tau_2 \dots d\tau_{M_s} $$
Vanishing of the residual vector **R** with a nonlinear solver means that the equilibrium steady state is reached. Usually, the iterative procedure takes the following form:

$$ \mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} - \mathbf{H}^{-1}\mathbf{R}(\mathbf{u}^{(k)}) \quad (10) $$

where **H** is a numerical estimation of the Jacobian matrix [9]. Vector $F_{nl}^u(\mathbf{u})$ in equation (8) is generally not an explicit function of the unknown vector **u** except for particular nonlinearities, such as cubic restoring forces for instance.

# 3 Application to simple cases

## 3.1 Algorithm and implementation

For the sake of generality, a program capable of taking into account any kind of nonlinearities in a systematic way has been developed. The cornerstone of the method is to accurately evaluate the nonlinear terms directly in the time domain. To this end, n-dimensional discrete direct (n-FFT) and inverse (n-iFFT) Fast Fourier Transforms are performed at each iteration of the nonlinear resolution in order to alternate between the frequency domain, where the nonlinear set of equations is solved, and the time domain, where the nonlinear terms are considered. Figure 1 illustrates the general algorithm. When the equations of motion are projected onto

<figure><img src="image_2.png" /><figcaption>Figure 1 – General implementation of the n-Harmonic Balance method</figcaption></figure>

the frequency domain, the notion of time has vanished. For a periodic signal, the hypertime domain collapses to one period. For an aperiodic behavior, it can be interesting to visualize true time going through the hypertime domain. Satisfying the relationship between true time

<figure><img src="image_3.png" /><figcaption>Figure 2 – Schematic of the path of time *t* through the hypertime domain [$\tau_1, \tau_2$]</figcaption></figure>
$t$ and pseudo-time variables $\tau_i$ yields in two dimensions:

$$ \left\{ \begin{array}{l} \tau_1 = \omega_1 t \\ \tau_2 = \omega_2 t \end{array} \right. \Rightarrow \tau_2 = \frac{\omega_2}{\omega_1} \tau_1 \qquad (11) $$

True time $t$ can be represented by the line $\tau_2 = \frac{\omega_2}{\omega_1} \tau_1$ in the $(\tau_1, \tau_2)$-plane. As it is periodic according to $\tau_1$ and $\tau_2$, it fulfils the light grey square $[0, 2\pi]^2$ depicted in figure 2. Indeed, the line segments numerated 2, 3 and 4 represent an interpretation of the mapping defined by: $t \in [0, +\infty[ \rightarrow (\tau_1, \tau_2) \in [0, 2\pi]^2$ such as equation (11) is verified. By construction, two points infinitely close to each other in $[0, 2\pi]^2$ may be completely different true instants in $[0, +\infty[$. This geometric interpretation can of course be generalized to a cubic hypertime domain or more.

## 3.2 Cubic nonlinearity

In order to check the implementation of the mentioned algorithm, a simply supported beam connected to a cubic nonlinear spring at $x = L/2$ as depicted in figure 3 (L is the length of the beam) is subjected to the following aperiodic external load $f_{ext} = 100(\cos(\pi t) + \cos(e^1 t))$. The displacement field is discretized by the usual Euler-Bernoulli functions. When the steady

<figure><img src="image_4.png" /><figcaption>Figure 3 – Schematic of the beam with a cubic restoring force</figcaption></figure>

state is reached, figure 4 confirms the accuracy of the HBM dealing with cubic nonlinearities. The maximum difference between time integration (namely the central finite difference scheme) results and 2-HBM results are of order $10^{-14}$ m. As explained in [10], only odd and odd-combination (odd sum of indexes in $\tau_1$ and $\tau_2$) harmonics appear in the solution. The residual error is obviously due to the truncation of the double Fourier series up to the number of harmonics $N_h$ kept. Trivially the higher the number of harmonics, the smaller the error. It has to be noted that a comparison of the CPU costs between the n-HBM and the time marching procedure is difficult: because of the aperiodic behavior of the motion and the mapping between true time $t$ and the hypertime domain, results given by the n-HBM are valid for $t \in [0, +\infty[$ that would require an infinite amount of time for the computation directly in the time domain. An analytical analysis is given below for comparison. For the sake of simplicity, the displacement is considered up to the first harmonic only ($\mathbf{u} = [a_1 \ a_2 \ b_1 \ b_2]^T$):

$$ v_{|x=L/2} = a_1 \cos \tau_1 + a_2 \cos \tau_2 + b_1 \sin \tau_1 + b_2 \sin \tau_2 \qquad (12) $$

Then, projecting the cubic nonlinear term onto the multi-frequency domain and performing the Galerkin procedure gives the respective components of $F_{nl}^u(\mathbf{u})$ (see table 1).

Reading table 2 shows the perfect agreement between the n-FFT technique and the analytical results for an explicit nonlinearity in the frequency domain. The precision of the alternating frequency/time domain approach is left to the reader's appreciation.
<figure><img src="image_5.png" /><figcaption>**Figure 4** – Displacement of the beam in $x = L/2$: time integration (—) and 2-HBM (- -) for $N_h = 7$</figcaption></figure>

<table><caption>**Table 1** – Analytical calculation of the cubic nonlinear restoring force expressed in the frequency domain for $N_h = 1$</caption><thead><tr><th>harmonic combinations</th><th>cubic force components</th></tr></thead><tbody><tr><td>cos τ<sub>1</sub></td><td>3k<sub>nl</sub>/4 (2a<sub>1</sub>b<sub>2</sub><sup>2</sup> + 2a<sub>2</sub><sup>2</sup>a<sub>1</sub> + a<sub>1</sub>b<sub>1</sub><sup>2</sup> + a<sub>1</sub><sup>3</sup>)</td></tr><tr><td>cos τ<sub>2</sub></td><td>3k<sub>nl</sub>/4 (a<sub>2</sub><sup>3</sup> + a<sub>2</sub>b<sub>2</sub><sup>2</sup> + 2a<sub>2</sub>a<sub>1</sub><sup>2</sup> + 2a<sub>2</sub>b<sub>1</sub><sup>2</sup>)</td></tr><tr><td>sin τ<sub>1</sub></td><td>3k<sub>nl</sub>/4 (2b<sub>2</sub><sup>2</sup>b<sub>1</sub> + a<sub>1</sub><sup>2</sup>b<sub>1</sub> + 2a<sub>2</sub><sup>2</sup>b<sub>1</sub> + b<sub>1</sub><sup>3</sup>)</td></tr><tr><td>sin τ<sub>2</sub></td><td>3k<sub>nl</sub>/4 (a<sub>2</sub><sup>2</sup>b<sub>2</sub> + b<sub>2</sub><sup>3</sup> + 2a<sub>1</sub><sup>2</sup>b<sub>2</sub> + 2b<sub>2</sub>b<sub>1</sub><sup>2</sup>)</td></tr></tbody></table>

<table><caption>**Table 2** – Comparison between 2-HBM results with (1) the cubic nonlinearity expanded directly in the multi-frequency domain and (2) using the n-FFT technique</caption><thead><tr><th rowspan="2">harmonics</th><th colspan="2">nonlinearity computed</th></tr><tr><th>directly in the frequency domain</th><th>using the n-FFT technique</th></tr></thead><tbody><tr><td>a<sub>1</sub></td><td>1.190871824484995 · 10<sup>-3</sup></td><td>1.190871824485002 · 10<sup>-3</sup></td></tr><tr><td>a<sub>2</sub></td><td>1.191004700829494 · 10<sup>-3</sup></td><td>1.191004700829497 · 10<sup>-3</sup></td></tr><tr><td>b<sub>1</sub></td><td>1.159230995878226 · 10<sup>-6</sup></td><td>1.159230995878233 · 10<sup>-6</sup></td></tr><tr><td>b<sub>2</sub></td><td>1.340061751905898 · 10<sup>-6</sup></td><td>1.340061751905895 · 10<sup>-6</sup></td></tr></tbody></table>

### 3.3 Contact constraint

The model of the previous section is considered now with a contact constraint (instead of a cubic spring) on $v_{|x=L/2}$ as depicted in figure 5. The new equation of motion becomes:

$$ \left\{ \begin{array}{l} \mathbf{M}\ddot{\mathbf{X}} + \mathbf{D}\dot{\mathbf{X}} + \mathbf{K}\mathbf{X} + \mathbf{F}_{\text{nl}}(\mathbf{X}, \dot{\mathbf{X}}) = \mathbf{F}_{\text{ext}} \\ v_{|x=L/2} \le v_{\text{max}}, \quad \forall t \ge t_0 \end{array} \right. \qquad (13) $$

The contact is treated as follows: the contact reactions are not explicitly calculated but if a penetration is detected in the time domain, the displacement $v_{|x=L/2}$ is corrected such as to be equal to $v_{max}$. This can be considered as an exact contact law where no penetration is allowed, equivalent to a dual formulation using Lagrange multipliers which is a straightforward extension. This approach is more accurate from a physical point of view than the one presented in [5] which only uses a penalty-based technique for taking contact constraints into account. This
<figure><img src="image_6.png" /><figcaption>**Figure 5** – Schematic of a beam constrained by a wall</figcaption></figure>

new displacement is then transformed back to the multi-frequency domain: this procedure is
depicted in figure 6. The 2-HBM results are compared to the time-integration results in figure 7.

<figure><img src="image_7.png" /><figcaption>**Figure 6** – Partial description of the 2-HBM algorithm dealing with a contact constraint. Displacement $v|_{x=L/2}$ is depicted. White harmonic combinations in subfigures (a) and (d) are redundant. For this specific case, the contact pressure is directly considered by correcting the displacement field in the hypertime domain as illustrated in subfigure (c)</figcaption></figure>

In order to get an accurate solution using the 2-HBM, a relatively high number of harmonics is needed. Unfortunately, the necessary number of spectral components that would guarantee an accurate solution is not known in advance. In fact, it is well known that all combinations of harmonics are more or less excited due the presence of the contact constraint as illustrated in figure 8 for which twelve harmonics for each frequency component (equivalent to a total of 313 non-redundant harmonic combinations) have been kept for the computation. Figure 7 shows that the solution from time integration and its counterpart from HBM are nearly indistinguishable. Further calculations would show that five harmonics for each frequency compo-
<figure><img src="image_8.png" /><figcaption>Figure 7 – Displacement $v_{|x|=L/2}$: time integration (—) and 2-HBM (- -) for $N_h = 12$</figcaption></figure>

nent represent a good compromise between CPU costs and accuracy. Indeed, the use of discret Fourier Transform techniques requires an estimation of $2N\lg_2 N$ computations if N stands for the number of evaluated points. Consequently, the number of computations greatly increases if a 2-dimensional (or more) domain is considered: this represents the main drawback of the proposed approach.

<figure><img src="image_9.png" /><figcaption>Figure 8 – Normalized absolute values of harmonic coefficients for $N_h = 12$</figcaption></figure>

Yet an effective incremental harmonic balance with two time scales for computing the response of the previous model with such a strong nonlinearity in time domain has been validated.

# 4 Rotor-stator modal interaction phenomenon

## 4.1 Basic background

Under certain conditions related to the rotational velocity of an aircraft engine, an energy exchange via direct structural contacts can lead to a wave speed coincidence between the $n$-nodal diameter modes of both structures. This very specific phenomenon is explained in de-
tails in [8]. It features a limit cycle associated to an aperiodic self-excited oscillating behavior of the entire structure even if any external load has disappeared [11]. Very simple geometrical and linear considerations state that a unique critical rotational velocity Ω can be dangerous in terms of amplitude of vibration:

$$ \frac{\omega_c}{n} = \Omega - \frac{\omega_{bd}}{n} \qquad (14) $$

where $\omega_{bd}$ and $\omega_c$ respectively stand for the bladed disk's and the casing's eigenfrequencies. In simple words, this condition stipulates that there is a rotational velocity for which the velocity in a stationary frame of a resonant backward wave travelling on the bladed disk ($\Omega - \frac{\omega_{bd}}{n}$) equals the velocity of a resonant forward travelling wave on the casing ($\frac{\omega_c}{n}$). An illustration of this phenomenon for $n = 3$ is given in figure 9. The black square ■ follows the forward wave propagating on the casing whereas the black blade illustrates the motion of the bladed disk. The aperiodic behavior is shown by the fact that after one round of the wave on the casing, the black blade is away from its initial position.

However, a numerical investigation using a finite element model including a contact law in conjunction with a time marching procedure showed that equation (14) is too restrictive and cannot precisely predict the interaction phenomenon [11]. Moreover, in order to overcome both the influence of initial conditions and the time-step size problem due to the contact detection, the frequency-domain method described in this paper can be helpful in giving a intrinsic description of the phenomenon.

<figure><img src="image_10.png" /><figcaption>Figure 9 – Snapshots of the bladed disk - casing model after three periods $T_c$</figcaption></figure>

## 4.2 Autonomous systems

Because the limit cycle explained above and which is of interest in this study occurs even when any external forcing has disappeared, it is solution of an autonomous set of equations. Tackling such a system of nonlinear equations, where the zero solution is trivial, relies on a modified version of the generalized harmonic balance described. In order to find a non-trivial aperiodic solution, $M_s$ incommensurate frequencies $\omega_m$ are included in the unknown vector. The well-posedness of the problem holds by prescribing $M_s$ harmonic coefficients among the u components. Because the system is autonomous, this has no influence on the almost periodic response except in the appearance of the phase shift. Equation (8) is then rewritten in the following manner:

$$ \mathbf{R}(\mathbf{u}, \omega) = \mathbf{H}(\omega)\mathbf{u} + \mathbf{F}_{nl}^{u}(\mathbf{u}, \omega) \qquad (15) $$

such has the unknown vector $\mathbf{u}$ is diminished of $M_s$ harmonic coefficients that are replaced by the $M_s$ incommensurate frequencies to be sought in order to find a non-zero solution and
arranged in ω. Convergence is obtained when the residual vector **R**(**u**, **ω**) is smaller than a given value.

# 5 Conclusions and prospects

The harmonic balance method has been adapted to solve systems with strong nonlinearities in a systematic way. Two simple cases, namely (1) a beam connected to a cubic nonlinear spring and (2) a beam constrained by a wall, showed the efficiency and accuracy of the *n*-dimensional HBM. Work are in progress to apply this method to a rotor-stator system undergoing an aperi- odic motion under very particular circumstances. It is believed that this technique will give a description of the motion close to the concept of nonlinear normal mode.

# References

[1] R.J. Gilmore and M.B. Steer. Nonlinear circuit analysis using the method of harmonic balance - a review of the art. part i. introductory concepts. *International Journal on Microwave and Millimeter Wave Computer Aided Engineering*, 1(1):22-37, 1991.

[2] S.L. Lau, Y.K. Cheung, and S. Y. Wu. Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems. *Journal of Applied Mechanics*, 50: 871-876, 1983.

[3] T. M. Cameron and J. H. Griffin. An alternating frequency/time domain method for calculating the steady state response of nonlinear dynamic systems. *Journal of Applied Mechanics*, 56(1):149-154, march 1989.

[4] A. Ushida and L.O. Chua. Frequency-domain analysis of nonlinear circuit driven by multi-tone signals. *IEEE Transactions on circuits and systems*, 31(9):766-779, 1984.

[5] Y.B. Kim and S.-K. Choi. A multiple harmonic balance method for the internal resonant vibration of a non-linear jeffcott rotor. *Journal of Sound and Vibration*, 208(5):745-761, 1997.

[6] R.R. Pušenjak and M. M. Oblak. Incremental harmonic balance method with multiple time variables for dynamical systems with cubic non-linearities. *International Journal for Numerical Methods in Engineering*, 59:255-292, 2004.

[7] M. Legrand. Modèles de prédiction de l'interaction rotor/stator dans un moteur d'avion. Thèse de doctorat, École Centrale de Nantes, 2005. oai:tel.archives-ouvertes.fr:tel-00011631_v2.

[8] P. Schmiechen. *Travelling Wave Speed Coincidence*. Ph.d. thesis, College of Science, Technology and Medicine, London, UK, 1997.

[9] R. Flechter. *Practical Methods of Optimization*. John Wiley, 1987.

[10] A.H. Nayfeh and B. Balachandran. *Applied Nonlinear Dynamics: analytical, computational and experimental methods*. Wiley Interscience, 1995.

[11] M. Legrand, B. Peseux, and C. Pierre. étude de l'interaction modale rotor/stator dans un moteur d'avion. In 14ᵉ *Colloque Vibrations, Chocs et Bruit*, Lyon, juin 2004.