# Une approche pour la modélisation et le contrôle des instabilités de combustion

Ioan D. Landau a,\* — Fethi Bouziani a — Robert. R. Bitmead b

a Laboratoire d'Automatique de Grenoble, ENSIEG BP 46,
38402 Saint-Martin d'Hères, France

b Department of Mechanical & Aerospace Engineering, University of California,
San Diego, La Jolla CA 92093-0411, USA

\* Corresponding author : Ioan-Dore.landau@inpg.fr

**RÉSUMÉ.** Un système de deux équations de Van der Pol généralisées et couplées est proposé comme un modèle de commande pour l'instabilité de combustion. Ce système est analysé en utilisant la méthode de Krylov-Bogoliubov. Les aspects de commande conduisant à l'extinction des oscillations sont examinés. Les résultats de l'analyse sont comparés avec des tests en simulation.

**ABSTRACT.** A set of two coupled generalized Van der Pol equations is proposed as a control model for combustion instabilities. The system is analyzed using the Krylov-Bogoliubov method. The control aspects related to quenching of the oscillations are examined. The analysis results are compared with simulation results.

**MOTS-CLÉS** : modélisation, instabilité de combustion, systèmes non-linéaires oscillants, méthode de Krylov-Bogoliubov.

**KEYWORDS** : modelling, combustion instabilities, nonlinear oscillating systems, Krylov-Bogoliubov method.
# 1. Prélude

Claude Lobry a apporté des contributions majeures à l'Automatique Non-linéaire et a contribué à l'essor de cette discipline en France. En plus de ses contributions scientifiques il a aussi transmis à ses collègues quelques message très importants :

– Un bon résultat de recherche sera reconnu par la communauté scientifique même si le texte est écrit en français.

– Il y a une "élégance cachée" dans un bon résultat de recherche. Il faut savoir le mettre en évidence.

# 2. Introduction

Le problème des instabilités de la combustion dans les turbines à gaz constitue l'objet de très nombreuses études. Il s'agit d'un phénomène complexe mais qui dans beaucoup de cas peut être expliqué par une contre réaction positive thermo-acoustique (Figure 1). Des programmes de recherche importants orientés sur cette thématique sont conduits dans un certain nombre de pays en collaboration avec les industriels (France, Etats-Unis, Grande-Bretagne). Les articles [1, 2, 3, 4, 5, 6, 7, 8] donnent une image de cette activité.

Le problème devient de plus en plus actuel car la réduction du rapport combustible/air, améliore le rendement et réduit la pollution. Malheureusement, la réduction du rapport combustible/air conduit à l'apparition des instabilités de combustion qui sont rédhibitoires. Plusieurs approches pour traiter ce problème ont été considérées. La première approche dite "passive" concerne la modification de la géométrie des turbines.

La seconde approche "active instability control" (AIC) est basée sur la mise en place d'une régulation en boucle fermée. Dans la plus part des cas, on traite l'instabilité comme une perturbation qu'on essaye de compenser. Ceci requiert la mise en place de capteurs et surtout d'actionneurs. Dans plusieurs applications il s'agit d'un haut parleur [9]. Cette approche de point de vue automatique est à rapprocher avec le contrôle actif des vibrations dans les systèmes mécaniques ou avec le contrôle de bruit. Cette approche ignore totalement le modèle du processus de combustion.

<figure><img src="image_1.png" /><figcaption>**Figure 1.** Schéma bloc décrivant l'instabilité thermo-acoustique</figcaption></figure>

Par ailleurs il a été constaté expérimentalement que dans certain cas une faible modulation en haute fréquence du rapport combustible/air peut conduire à la disparition des instabilités de combustion (cette commande s'apparente à une entrée multiplicative dans
le système représenté dans Figure 1). La mise en œuvre systématique repose néanmoins sur la disponibilité d'un modèle du mécanisme d'instabilité et sur la compréhension des conditions assurant l'extinction des oscillations.

Cette approche, dans la mesure où un support théorique peut être apporté présente d'une part des avantages pratiques évidents en terme de simplicité de mise en œuvre et d'autre part en terme de rendement énergétique (pas d'apport d'énergie extérieure).

Plusieurs aspects importants sont à étudier :

1) Un modèle macroscopique du mécanisme de l'instabilité ;

2) Une compréhension théorique approfondie des comportements oscillatoires du modèle ;

3) Etude théorique des conditions d'extinction des oscillations par application d'une commande multiplicative.

Des progrès important ont été fait sur le point 1 dans [10, 11]. Un modèle assez représentatif existe. A titre de référence, nous avons considéré le modèle de Dunstan et Bitmead [12, 13] qui est une extension du modèle Peracchio et Proscia [14] et qui a fait l'objet de nombreuses vérifications expérimentales. Ce modèle est représenté dans la Figure 2, où $p_t$ est la perturbation de pression en aval agissant de la chambre de combustion, $q_t$ est la quantité de chaleur dégagée par la flamme, $\varphi[\cdot]$ est une fonction non linéaire statique, $\tau$ est un temps de retard du transport de l'injecteur à la surface de la flamme, $\omega$ est le mode fondamental acoustique de l'instabilité de combustion, $\zeta$ est le coefficient d'amortissement associé et N et M sont les gains des résonateurs de deuxième ordre avec les pulsations naturelles égale à $\omega$ et $3\omega$, respectivement. Le modèle résulte de l'interaction entre l'acoustique de la chambre de combustion et l'emplacement de la surface de la flamme.

La caractéristique statique non-linéaires $\varphi[\cdot]$ a été identifiée sur une turbine réelle et on a une image assez claire de ses caractéristiques [12].

<figure><img src="image_2.png" /><figcaption>**Figure 2. Modèle de l'instabilité de combustion**</figcaption></figure>

Le modèle est caractérisé par une structure à contre-réaction intrinsèque résultant de la contre réaction thermo-acoustique (voir Figure 2). La voie directe est caractérisée par la présence de deux résonateurs et le retour non-linéaire crée dans certaines situations un amortissement négatif conduisant à l'apparition des oscillations.
En faisant abstraction de l'existence de deux résonateurs et du retard, le schéma sug-
gère une analogie avec les oscillateurs électroniques. En effet, si en plus on examine la
caractéristique non-linéaire celle-ci peut être approchée dans la zone normale d'opération
par

$$
\varphi(y) = a + by - \frac{c}{3}y^3. \tag*{[1]}
$$

où a, b et c sont des constantes.

Ceci suggère un rapprochement avec l'équation de Van der Pol largement utilisée
pour l'étude des oscillateurs électroniques et en général avec les techniques d'études des
systèmes oscillants. Il faut néanmoins ne pas perdre de vue qu'à l'opposé des oscillateurs
électroniques l'objectif et de faire disparaître les oscillations et non pas de les amorcer.

En effet l'équation de Van der Pol a la forme

$$
\frac{d^2 y}{dt^2} + \omega^2 y = \epsilon (1 - y^2) \frac{dy}{dt}, \quad [2]
$$

où $\epsilon$ est paramètre faible, $\omega$ est la pulsation naturelle et dont le terme de droite correspond
à la dérivée de la fonction non-linéaire

$$
\varphi_v(y) = \varphi_{v0} + y - \frac{y^3}{3}, \quad [3]
$$

où $\varphi_{v0}$ est une constante arbitraire.

Tenant compte des caractéristiques de la non-linéarité (1) effective, une généralisation
de l'équation de Van der Pol doit être considérée (Van der Pol généralisée, VDPG)

$$
\frac{d^2 y}{dt^2} + \omega^2 y = \frac{d}{dt} \left( \varphi_{v0} + \varphi_{v1} y - \frac{\varphi_{v3}}{3} y^3 \right), \quad [4]
$$

où $\varphi_{v1}$ et $\varphi_{v3}$ sont des constantes arbitraires positive, et qui dans le cas de deux résonateurs conduit au schéma représenté dans la Figure 3 et au système d'équation suivant

$$
\left\{
\begin{array}{l}
\frac{d^2 x_1}{dt^2} + \omega_1^2 x_1 = \epsilon \frac{d}{dt} \left( \varphi_{v0} + \varphi_{v1} y - \frac{\varphi_{v3}}{3} y^3 \right), \\
\frac{d^2 x_2}{dt^2} + \omega_2^2 x_2 = \epsilon \frac{d}{dt} \left( \varphi_{v0} + \varphi_{v1} y - \frac{\varphi_{v3}}{3} y^3 \right), \\
y = x_1 + x_2.
\end{array}
\right.
\qquad [5]
$$

Ce modèle constitue sans aucun doute dans un premier temps, une base pour l'étude des
phénomènes d'oscillations (qui serait à confronter avec les données réelles) et dans un
deuxième temps il devrait permettre l'étude effective de l'effet d'extinction, et donc le
dimensionnement de la commande (amplitude et caractéristiques fréquentielles).

Une fois le modèle adopté, il s'agit de trouver l'approche mathématique permettant
d'étudier analytiquement ce modèle. Les équations régissant ce système, font que ce der-
nier appartient à la classe des systèmes qui sont décrits par des équations de la forme

$$
\frac{d^2 x}{dt^2} + \omega^2 x = \epsilon f \left( x, \frac{dx}{dt} \right), \quad [6]
$$

où $\epsilon$ est une quantité positive de valeurs faible et $f$ peut être une série entière de $\epsilon$ dont les
coefficients sont polynomial en $x$ et $\frac{dx}{dt}$.
<figure><img src="image_3.png" /><figcaption>**Figure 3. Modèle de l'instabilité de combustion basé sur des équations de Van der Pol généralisées et couplées**</figcaption></figure>

Puisqu'en général on ne peut pas trouver les solutions exactes pour ce type d'équation, les procédures d'approximation pour l'analyse de ce type d'équations doivent être considérées.

La méthode de Krylov-Bogoliubov (K-B) [15, 16, 17, 18, 19, 20] est sans doute une des procédures les plus efficaces pour analyser les systèmes oscillants régis par des équations de la forme (6). En bref la méthode K-B cherche des solutions de la forme

$$x(t) = a(t) \cos \psi(t), \qquad [7]$$

où $a$ est l'amplitude de l'oscillation fondamentale variant dans le temps et $\psi$ est la phase totale instantanée. Ils obéissent dans le cas d'un seul résonateur au système d'équations

$$\begin{cases} \frac{da}{dt} = -\frac{1}{2\omega\pi} \int_{0}^{2\pi} f(a \cos \psi, -a\omega \sin \psi) \sin \psi d\psi, \\ \frac{d\psi}{dt} = \omega - \frac{1}{2\pi\omega a} \int_{0}^{2\pi} f(a \cos \psi, -a\omega \sin \psi) \cos \psi d\psi. \end{cases}$$

Notons que $\psi$ peut être écrite sous la forme

$$\psi(t) = \omega t + \theta(t),$$

où $\theta$ est la phase instantanée

C'est cette approche qui sera employée pour l'analyse du comportement du modèle d'instabilité de combustion. Du point de vue de l'analyse, le modèle d'instabilité de combustion présente un certain nombre de difficultés parmi lesquelles nous mentionnons :

- La présence de deux résonateurs couplés.

- La dynamique compliquée dans le retour de boucle qui est due à la présence en cascade d'une dérivée et d'un retard.

Le travail présent se concentrera sur l'analyse de l'effet de la structure de deux résonateurs couplés, sans que nécessairement il existe une relation entre leur harmoniques.
### 3. Première approximation K-B pour les systèmes multi-résonateurs autonomes

Considérons un système avec $n$ résonateurs qui est décrit par des équations de la forme

$$ \frac{d^2 x_k}{dt^2} + \omega_k^2 x_k = \epsilon f_k \left( x, \frac{dx}{dt} \right), \quad (k = 1, 2, \dots, n) \qquad [8] $$

où $x = \{x_1, \dots, x_n\}$, $\frac{dx}{dt} = \{\frac{dx_1}{dt}, \dots, \frac{dx_n}{dt}\}$ et $\epsilon$ est un paramètre faible.

Pour récapituler (pour plus de détails voir le Chapter 2 de [18]), pour le résonateur $j$, la première approximation K-B propose la solution

$$ x_j = a_j \cos(\psi_j), \qquad [9] $$

où $\psi_j = \omega_j t + \theta_j$, $a_j$ et $\theta_j$ sont des fonctions lentement variant dans le temps obéissent aux équations

$$ \left\{ \begin{array}{l} \frac{da_j}{dt} = -\frac{\epsilon}{2\omega_j} H_{jj}(a_1, \dots, a_n, \theta_1, \dots, \theta_n), \\ \frac{d\theta_j}{dt} = -\frac{\epsilon}{2\omega_j a_j} G_{jj}(a_1, \dots, a_n, \theta_1, \dots, \theta_n). \end{array} \right. \qquad [10] $$

avec $H_{jj}$ et $G_{jj}$ sont obtenues de la fonction $f_j(x, \frac{dx}{dt})$ en injectant

$$ \left\{ \begin{array}{l} x_k = a_k \cos(\omega_k t + \theta_k), \\ \frac{dx_k}{dt} = -a_k \omega_k \sin(\omega_k t + \theta_k), \end{array} \right. \quad (k = 1, 2, \dots, n) \qquad [11] $$

et en la mettant sous la forme

$$ \begin{aligned} f_j & (a_1 \cos(\omega_1 t + \theta_1), \dots, a_n \cos(\omega_n t + \theta_n), \\ & -a_1 \omega_1 \sin(\omega_1 t + \theta_1), \dots, -a_n \omega_n \sin(\omega_n t + \theta_n)) \\ = & H_{jj} \sin(\omega_j t + \theta_j) + G_{jj} \cos(\omega_j t + \theta_j) \\ & + \sum_{\substack{\ell=1 \\ \ell \neq j}}^{r} (H_{\ell j} \sin(\omega_{\ell} t + \theta_{\ell}) + G_{\ell j} \cos(\omega_{\ell} t + \theta_{\ell})), \end{aligned} \qquad [12] $$

où $\omega_\ell$ et $\theta_\ell$ sont des combinaisons linéaires de $\omega_1, \dots, \omega_n$ et $\theta_1, \dots, \theta_n$, respectivement, et $r$ est le nombre combinaisons linéaires possible de $\omega_1, \dots, \omega_n$ différent de $\omega_j$. En outre pour $x_j$, les coefficients du terme fondamental dans (12) sont utilisés et toutes les autres termes sont éliminés.

### 4. Application de l'approximation K-B pour deux équations de Van der Pol généralisées

Considérons le système d'équations (5) et la forme (8), dans ce cas

$$ \begin{aligned} f_1 &= f_2 = f(x_1, x_2, \frac{dx_1}{dt}, \frac{dx_2}{dt}) \\ &= \varphi_{v0} + \varphi_{v1} \left(1 - \frac{\varphi_{v3}}{\varphi_{v1}} (x_1 + x_2)^2\right) \left(\frac{dx_1}{dt} + \frac{dx_2}{dt}\right). \end{aligned} \qquad [13] $$

Introduisons
$$
\begin{cases}
x_i = a_i \cos(\omega_i t + \theta_i), \\
\frac{dx_i}{dt} = -a_i \omega_i \sin(\omega_i t + \theta_i),
\end{cases}
\quad (i = 1, 2)
$$

dans (13), on obtient

$$
\begin{align*}
& f(a_1 \cos(\omega_1 t + \theta_1), a_2 \cos(\omega_2 t + \theta_2), \\
& \qquad -a_1 \omega_1 \sin(\omega_1 t + \theta_1), -a_2 \omega_2 \sin(\omega_2 t + \theta_2)) \\
& = -\varphi_{v1} \left(1 - \frac{\varphi_{v3}}{\varphi_{v1}} (a_1 \cos(\omega_1 t + \theta_1) + a_2 \cos(\omega_2 t + \theta_2))^2\right) \\
& \quad \times (a_1 \omega_1 \sin(\omega_1 t + \theta_1) + a_2 \omega_2 \sin(\omega_2 t + \theta_2)) + \varphi_{v0}. \tag{14}
\end{align*}
$$

Pour approximer la solution de (5), il est nécessaire de mettre (14) sous la forme (12).
Dans [11], on donne les détails du calcul qui nous permet d’obtenir l’expression

$$
\begin{align*}
& f(a_1 \cos(\omega_1 t + \theta_1), a_2 \cos(\omega_2 t + \theta_2), \\
& \quad -a_1 \omega_1 \sin(\omega_1 t + \theta_1), -a_2 \omega_2 \sin(\omega_2 t + \theta_2)) \\
& = \varphi_{v0} + \varphi_{v1} \left\{ \begin{aligned}
& -\omega_1 a_1 \left(1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left(\frac{a_1^2}{4} + \frac{a_2^2}{2}\right)\right) \sin(\omega_1 t + \theta_1) \\
& -\omega_2 a_2 \left(1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left(\frac{a_2^2}{4} + \frac{a_1^2}{2}\right)\right) \sin(\omega_2 t + \theta_2) \\
& + \frac{\varphi_{v3}}{\varphi_{v1}} \left[ \omega_1 \frac{a_1^3}{4} \sin(3(\omega_1 t + \theta_1)) + \omega_2 \frac{a_2^3}{4} \sin(3(\omega_2 t + \theta_2)) \right. \\
& \qquad \left. + (2\omega_1 + \omega_2) \frac{a_1^2 a_2}{2} \sin((2\omega_1 + \omega_2)t + 2\theta_1 + \theta_2) \right. \\
& \qquad \left. + (\omega_1 + 2\omega_2) \frac{a_1 a_2^2}{2} \sin((\omega_1 + 2\omega_2)t + \theta_1 + 2\theta_2) \right. \\
& \qquad \left. + (2\omega_1 - \omega_2) \frac{a_1^2 a_2}{4} \sin((2\omega_1 - \omega_2)t + 2\theta_1 - \theta_2) \right. \\
& \qquad \left. + (2\omega_2 - \omega_1) \frac{a_2^2 a_1}{4} \sin((2\omega_2 - \omega_1)t + 2\theta_2 - \theta_1) \right]
\end{aligned} \right\}. \tag{15}
\end{align*}
$$

En examinant cette fonction on peut constater l'existence de l'ensemble de fréquences

$$
W = \{ \omega_1, \omega_2, 3\omega_1, 3\omega_2, 2\omega_1 + \omega_2, \omega_1 + 2\omega_2, \\
\qquad 2\omega_1 - \omega_2, 2\omega_2 - \omega_1 \}. \tag{16}
$$

Cet ensemble est très important pour trouver les différents régimes possibles du système,
c-a-d pour x₁ (respectivement x₂), les termes restants de (15) après l'application de l'ap-
proximation K-B serons seulement les termes avec la fréquence ω de W tell que ω ≈ ω₁
(respectivement ω₂). En conséquence, on a la classification suivante :

1) $\omega_1 \not\approx \{\omega_2, 3\omega_2, \frac{\omega_2}{3}\}$-deux générateurs avec extinction compétitive

2) $\omega_1 \approx \omega_2$-synchronisation mutuelle avec des fréquences proches

3) $\omega_1 \approx 3\omega_2$ (respectivement $\omega_2 \approx 3\omega_1$)-synchronisation mutuelle avec des fréquences multiples

4.1. Deux générateurs avec extinction compétitive

Considérons le cas où les fréquences $\omega_1$ et $\omega_2$ respectent la condition 1 cité ci-dessus.
Dans ce cas, ils n'existe aucun effet d'interconnexion entre les deux fréquences et l'ap-
proximation K-B utilise seulement les termes de l'oscillation fondamental de

$$
\begin{align*}
& f(a_1 \cos(\omega_1 t + \theta_1), a_2 \cos(\omega_2 t + \theta_2), \\
& \quad -a_1 \omega_1 \sin(\omega_1 t + \theta_1), -a_2 \omega_2 \sin(\omega_2 t + \theta_2)).
\end{align*}
$$
En conséquence les solutions approximées de (5) sont (pour les détails voir [11])

$$x_i = a_i \cos(\omega_i t + \theta_i), \quad (i = 1, 2) \qquad [17]$$

avec

$$\left\{ \begin{aligned} \frac{da_1}{dt} &= \epsilon \varphi_{v1} \frac{a_1}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_1^2}{4} + \frac{a_2^2}{2} \right) \right), \\ \frac{da_2}{dt} &= \epsilon \varphi_{v1} \frac{a_2}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_2^2}{4} + \frac{a_1^2}{2} \right) \right), \\ \frac{d\theta_1}{dt} &= 0, \\ \frac{d\theta_2}{dt} &= 0. \end{aligned} \right. \qquad [18]$$

Trouvons les solutions d'équilibre de (18). Dans ce cas (18) possède quatre solutions d'équilibre :

$$a_1 = 0 \text{ et } a_2 = 0, \qquad [19]$$

$$a_1 = \frac{2}{\sqrt{3}} \sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}} \text{ et } a_2 = \frac{2}{\sqrt{3}} \sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}, \qquad [20]$$

$$a_1 = 2 \sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}} \text{ et } a_2 = 0, \qquad [21]$$

$$a_1 = 0 \text{ et } a_2 = 2 \sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}. \qquad [22]$$

Les deux premières solutions (19) et (20) sont instables, et les deux dernières solutions (21) et (22) sont stables. Donc les amplitudes de $x_1$ et $x_2$ convergent à un des deux états stationnaires possibles (21) et (22). Selon la condition initiale, un des deux générateurs est excité, tandis que les oscillations de l'autre générateur sont entièrement éteintes. Une telle extinction des oscillations d'un des générateurs, causée par la présence d'un très fort couplage non linéaire entre les équations, est connue sous le nom *extinction compétitive*.

Il a été observé que si $a_1(0) > a_2(0)$, $x_1$ oscille et les oscillations de $x_2$ sont entièrement éteintes, et l'effet inverse se produit quant $a_1(0) < a_2(0)$. La Figure 4 présente un test de simulation avec $\omega_1 = 2\pi \times 210$, $\omega_2 = 2\pi \times 740$, $\epsilon = 1$, $\varphi_{v1} = 1.78 \times 10^2$, $\varphi_{v3} = 1.24 \times 10^7$ (valeurs approchant les valeurs réelles, voir [12]), $a_1(0) = 4 \times 10^{-3}$ et $a_2(0) = 2 \times 10^{-3}$, la partie supérieure montre les sorties simulées de (5) et la partie inférieure montre les sorties approximées par (17) et (18). Aussi, quant $a_1(0) = a_2(0)$ (mais non nuls), il a été observé que :

- Dans (18), les amplitudes $a_1$ et $a_2$ convergent à $\frac{2}{\sqrt{3}}\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}$ et $\frac{2}{\sqrt{3}}\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}$ respectivement, qui correspond à la solution équilibre instable (20).

- Dans (5), les amplitudes $a_1$ et $a_2$ convergent $\frac{2}{\sqrt{3}}\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}$ et $\frac{2}{\sqrt{3}}\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}}$ respectivement, et reste temporairement, mais après un lent moment (si on le compare à la dynamique de convergence) ces amplitudes vont converger nécessairement à une des deux solutions d'équilibres (21) et (22).

Cela implique que dans certaines conditions, les deux fréquences peuvent coexister pour un long moment avant d'entrer dans le régime d'extinction compétitive. Pour illustrer ce phénomène, la Figure 5 présente un test de simulation avec $\omega_1 = \pi$, $\omega_2 = 3.5\omega_1 = 3.5\pi$, $\epsilon = 0.1$, $\varphi_{v1} = \varphi_{v3} = 1$, $a_1(0) = 1$ et $a_2(0) = 1$, la partie supérieure est la sortie $x_1$ de (5) et la partie inférieure est la sortie $x_2$ de (5).
<figure><img src="image_4.png" /><figcaption>**Figure 4.** (A) $x_1$ simulé de (5), (B) $x_2$ simulé de (5), (C) $x_1$ approximé par (17) et (18), (D) $x_2$ approximé par (17) et (18)</figcaption></figure>

## 4.2. Synchronisation mutuelle avec des fréquences proches

Considérons le cas où les fréquences $\omega_1$ et $\omega_2$ sont proches. Pour $x_1$ (respectivement $x_2$), l'application de l'approximation K-B implique la conservation de tous les coefficients des termes sinusoïdaux dans $f(a_1 \cos(\omega_1 t + \theta_1), a_2 \cos(\omega_2 t + \theta_2), -a_1\omega_1 \sin(\omega_1 t + \theta_1), -a_2\omega_2 \sin(\omega_2 t + \theta_2))$ avec une fréquence proche de $\omega_1$ (respectivement $\omega_2$) et l'élimination de tous les autres termes. En conséquence les solutions approximées de (5) sont (voir [11] pour les détails)

$$x_i = a_i \cos(\omega_i t + \theta_i), \quad (i = 1, 2) \qquad [23]$$
avec $a_1, a_2, \theta_1$ et $\theta_2$ sont gouverné par

$$
\left\{
\begin{aligned}
\frac{da_1}{dt} &= \epsilon\varphi_{v1} \left\{ \frac{a_1}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_1^2}{4} + \frac{a_2^2}{2} \right) \right) \right. \\
&\qquad \left. + \left[ \frac{a_2\omega_2}{2\omega_1} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_1^2+a_2^2}{4} \right) - \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_2a_1^2}{4} \right] \cos(\Delta\psi) \right. \\
&\qquad \left. + (\omega_1 - 2\omega_2) \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_1a_2^2}{8\omega_1} \right) \cos(2\Delta\psi) \right\}, \\[2em]
\frac{da_2}{dt} &= \epsilon\varphi_{v1} \left\{ \frac{a_2}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_2^2}{4} + \frac{a_1^2}{2} \right) \right) \right. \\
&\qquad \left. + \left[ \frac{a_1\omega_1}{2\omega_2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_1^2+a_2^2}{4} \right) - \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_1a_2^2}{4} \right] \cos(\Delta\psi) \right. \\
&\qquad \left. + (\omega_2 - 2\omega_1) \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_2a_1^2}{8\omega_1} \right) \cos(2\Delta\psi) \right\}, \\[2em]
\frac{d\theta_1}{dt} &= -\epsilon\varphi_{v1} \left\{ (\omega_1 - 2\omega_2) \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_2^2}{8\omega_1} \sin(2\Delta\psi) \right. \\
&\qquad \left. + \left[ \frac{a_2\omega_2}{2a_1\omega_1} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_2^2}{4} + \frac{a_1^2}{2} \right) \right) \right. \right. \\
&\qquad \left. \left. + (2\omega_1 - \omega_2) \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_2a_1}{8\omega_1} \right] \sin(\Delta\psi) \right\}, \\[2em]
\frac{d\theta_2}{dt} &= \epsilon\varphi_{v1} \left\{ \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_1^2}{8\omega_2} (\omega_2 - 2\omega_1) \sin(2\Delta\psi) \right. \\
&\qquad \left. + \left[ \frac{a_1\omega_1}{2a_2\omega_2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_1^2}{4} + \frac{a_2^2}{2} \right) \right) \right. \right. \\
&\qquad \left. \left. + \frac{\varphi_{v3}}{\varphi_{v1}} \frac{a_2a_1}{8\omega_2} (2\omega_2 - \omega_1) \right] \sin(\Delta\psi) \right\}.
\end{aligned}
\right.
\quad [24]
$$

où $\Delta\psi = \psi_1 - \psi_2 = (\omega_1 - \omega_2)t + \theta_1 - \theta_2$.

<figure><img src="image_5.png" /><figcaption>**Figure 5.** (A) $x_1$ simulé de (5), (B) $x_2$ simulé de (5)</figcaption></figure>

Ce résultat est très important, car en parallèle de l'équation différentielle (5), il est possible de suivre l'évolution de l'amplitude et la phase de la sortie de (24) et de comparer les deux signaux qui sont mesurables en pratique.

L'intégration et l'étude des solutions stationnaires de (24) sont très difficiles. Cependant, pour trouver les solutions stationnaires quant $\omega_1 = \omega_2 = \omega$, on peut adopter les étapes suivantes.
Utilisons $y = x_1 + x_2$, si on additionne les deux équations de (5) on obtient

$$ 
\begin{aligned}
& \frac{d^2(x_1+x_2)}{dt^2} + \omega^2(x_1+x_2) \\
&= \epsilon\varphi_{v1}\left(1 - \frac{\varphi_{v3}}{\varphi_{v1}}(x_1+x_2)^2\right) \frac{d(x_1+x_2)}{dt} \\
&\Rightarrow \frac{d^2y}{dt^2} + \omega^2y = \epsilon\varphi_{v1}\left(1 - \frac{\varphi_{v3}}{\varphi_{v1}}y^2\right) \frac{dy}{dt}.
\end{aligned}
\quad [25] $$

On voit que (25) correspond à une équation de Van der Pol généralisée. Or on sait bien que pour une seule équation de Van der Pol généralisée, l'approximation K-B donne une solution stationnaire [19, 15, 16, 17, 18, 20]

$$ y = 2\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}} \cos(\omega t + \theta). \quad [26] $$

$$ \Rightarrow x_1 + x_2 = 2\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}} \cos(\omega t + \theta), \quad [27] $$

où $\theta$ est le déphasage entre $x_1$ et $x_2$. Et donc $a_1, a_2$ et $\theta_1 - \theta_2$ doivent satisfaire

$$ a_1^2 + a_2^2 + 2a_1a_2 \cos(\theta_1 - \theta_2) = 4 \frac{\varphi_{v1}}{\varphi_{v3}} \quad [28] $$

Notons que pour ce cas, il existe une infinité de points d'équilibre et que la convergence de l'amplitude et de la phase dépend essentiellement de l'état initial de $x_1$ et $x_2$. Par conséquence, pour avoir le même résultat (5) et (24), il faut initialiser (24) avec les bonnes valeurs d'amplitude et de phase.

<figure><img src="image_6.png" /><figcaption>**Figure 6.** (A) $x_1$ simulé de (5), (B) $x_2$ simulé de (5), (C) $x_1$ approximé par (23) et (24), (D) $x_2$ approximé par (23) et (24)</figcaption></figure>

La Figure 6 montre un test de simulation avec $\omega_1 = \omega_2 = 2\pi \times 210$, $\varphi_{v1} = 1.78 \times 10^2$, $\varphi_{v3} = 1.24 \times 10^7$, $a_1(0) = 2 \times 10^{-3}$, $a_2(0) = 10^{-3}$, $\epsilon = 1$ et $\theta_1(0) = \theta_1(0) = 0$, la partie supérieure montre les sorties simulées de (5) et la partie inférieure montre les sorties approximées par (23) et (24).
4.3. Synchronisation mutuelle avec des fréquences multiples

Considérons le cas où la fréquence $\omega_1$ est proche de $3\omega_2$. Dans ce cas les termes avec les fréquences $\omega_1$ et $3\omega_2$ sont utilisés pour l'approximation de $x_1$, et les termes avec les fréquences $\omega_2$ et $(2\omega_2 - \omega_1)$ sont utilisés l'approximation de $x_2$. Donc, on trouve (voir [11] pour les détails)

$$x_i = a_i \cos(\omega_i t + \theta_i), \quad (i = 1, 2)$$

avec

$$
\left\{
\begin{aligned}
\frac{da_1}{dt} &= \epsilon \varphi_{v1} \left[ \frac{a_1}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_1^2}{4} + \frac{a_2^2}{2} \right) \right) \right. \\
&\qquad \left. - \frac{\varphi_{v3}}{\varphi_{v1}} \frac{\omega_2 a_2^3}{8\omega_1} \cos(\Delta\psi) \right], \\
\frac{da_2}{dt} &= \epsilon \varphi_{v1} \left[ \frac{a_2}{2} \left( 1 - \frac{\varphi_{v3}}{\varphi_{v1}} \left( \frac{a_2^2}{4} + \frac{a_1^2}{2} \right) \right) \right. \\
&\qquad \left. - \frac{a_1 a_2}{8\omega_2} \frac{\varphi_{v3}}{\varphi_{v1}} (\omega_1 - 2\omega_2) \cos(\Delta\psi) \right], \\
\frac{d\theta_1}{dt} &= \epsilon \varphi_{v3} \frac{\omega_2 a_2^3}{8\omega_1 a_1} \sin(\Delta\psi), \\
\frac{d\theta_2}{dt} &= -\epsilon \varphi_{v3} \frac{a_1 a_2}{8\omega_2} (\omega_1 - 2\omega_2) \sin(\Delta\psi),
\end{aligned}
\right.
$$

$$\text{où } \Delta\psi = \psi_1 - 3\psi_2 = (\omega_1 - 3\omega_2)t + \theta_1 - 3\theta_2.$$

Pour (29) il existe deux points d'équilibre stables, un point peut être calculé analytiquement de (29)

$$a_1 = 2\sqrt{\frac{\varphi_{v1}}{\varphi_{v3}}} \quad \text{et} \quad a_2 = 0. \qquad [29]$$

Le deuxième point qui concerne le phénomène de synchronisation a besoin d'une solution numérique de (29). Par exemple pour le cas où $\omega_1 = 3\omega_2 = 6\pi \times 210$, $\epsilon = 1$, $\varphi_{v1} = 1.78 \times 10^2$ et $\varphi_{v1} = 1.24 \times 10^7$, le deuxième point d'équilibre correspond à

$$a_1 = 2.25 \times 10^{-3}, \quad a_2 = 8.09 \times 10^{-3} \quad \text{et} \quad \Delta\psi = \pi. \qquad [30]$$

De cela, on peut voir que si $\omega_1$ est près de $3\omega_2$ (respectivement $\omega_2/3$), il est possible d'avoir deux phénomènes selon l'état initial. Dans le premier phénomène, le générateur avec la fréquence $\omega_1$ est excité et l'autre générateur avec la fréquence $\omega_2$ est éteint. Dans le deuxième phénomène, on a le régime de synchronisation. On entend par le régime de synchronisation, un régime où la fréquence des oscillations du deuxième générateur qui est égale à $\omega_2 + \dot{\theta}_2$, est exactement un tiers de la fréquence des oscillations du premier générateur qui est égale à $\omega_1 + \dot{\theta}_1$.

La Figure 7 présente un test de simulation avec $\omega_1 = 3\omega_2 = 6\pi \times 210$, $\epsilon = 1$, $\varphi_{v1} = 1.78 \times 10^2$, $\varphi_{v3} = 1.24 \times 10^7$, $a_1(0) = 10^{-3}$, $a_2(0) = 6 \times 10^{-3}$ et $\theta_1(0) = \theta_1(0) = 0$, la partie supérieure montre les sorties simulées de (5) et la partie inférieure montre les sorties approximées par (29) et (29).
<figure><img src="image_7.png" /><figcaption>**Figure 7.** (A) $x_1$ simulé de (5), (B) $x_2$ simulé de (5), (C) $x_1$ approximé par (29) et (29), (D) $x_2$ approximé par (29) et (29)</figcaption></figure>

**5. Extinction des oscillations**

Du point de vue pratique, il est possible de moduler le flux du combustible entrant dans la chambre de combustion. Cela se traduit par une excitation multiplicative $r_t$ agissant sur le retour de la boucle du modèle de l'instabilité de combustion. Si on considère le modèle de Dunstan/Bitmead (Figure 2) cela se traduit par la modification montrée dans la Figure 8.

<figure><img src="image_8.png" /><figcaption>**Figure 8.** Effet de la modulation sur modèle de Dunstan et Bitmead</figcaption></figure>

Une étude sur l'effet de cette excitation a été réalisée en considérant un modèle simplifié avec un seul résonateur (c-à-d un seul résonateur de VDPG). Le but de cette étude était de voir si on peut arriver à un phénomène d'extinction des oscillations dans le système de combustion en utilisant une excitation $r_t = B \cos(w_r t)$ en hautes fréquences.

En effet, l'étude a démontré que l'extinction des oscillations est possible sous certaines conditions. Un exemple de simulation sur l'extinction des oscillations est montré dans la
Figure 9 avec $w = 2\pi \times 210$, $w_r = 6\pi \times 2100$, $\varphi_{v0} = 0.45$, $\varphi_{v1} = 1.78 \times 10^2$, $\varphi_{v3} = 1.24 \times 10^7$ et un filtre passe bas de premier ordre. L'étude en question fera l'objet d'une publication ultérieure.

<figure><img src="image_9.png" /><figcaption>Figure 9. Exemple de simulation sur l'extinction des oscillations</figcaption></figure>

# 6. Conclusion

Le but de cet article a été de prouver qu'il est possible d'aller plus loin dans l'analyse du modèle d'instabilité de combustion proposé dans [12] et d'envisager le problème de commande des instabilités. Les deux équations Van der Pol généralisées et couplées considérées dans cet article peuvent être un choix efficace pour approcher le modèle d'instabilité de combustion. La méthode d'analyse est basée sur l'utilisation de l'approche de Krylov-Bogoliubov pour les systèmes oscillants. En effet, cette approche nous a permis de surmonter une des difficultés liées au modèle d'instabilité de combustion, c-à-d la présence de deux résonateurs couplés. Les essais de simulation ont illustré la précision de l'approximation de Krylov-Bogoliubov.

# 7. Bibliographie

[1] V. Faivre et T. Poinsot. Experimental and numerical investigations of jet active control for combustion applications. *Journal of Turbulence*, 5(025), Aug 2004.

[2] R. Blonbou et A. Laverdant. Control of combustion instabilities on a rijke tube by a neural network. TP 49, ONERA, 29, avenue de la Division Leclerc, BP 72, 92322 Chltillon Cedex, France, May 2000.
[3] C.A. Jacobson, A. Khibnik, A. Banaszuk, J.M. Cohen, et W.P. Proscia. Active control of com-
bustion instabilities in gas turbine engines for low emissions. part I : Physics-based and experi-
mentally identified models of combustion instability. *Proc. AVT Symposium on Active Control*
*Technology, Braunschweig*, pages 30–1–30–11, 2000.

[4] K. McManus, F. Han, W. Dunstan, C. Barbu, et M. Shah. Modeling and control of combustion dynamics in industrial gas turbines. *Proc. ASME Turbo-Expo*, pages 567–575, 2004.

[5] E. Haile M. Mettenleiter et S. Candel. Adaptive control of aeroacoustic instabilities. *J. Sound & Vibration*, 230(4):761–789, 2000.

[6] A.P. Dowling. The calculation of thermoacoustic oscillations. *J Sound & Vibration*, 180(4):557–581, 1995.

[7] S. Candel. Combustion dynamics and control : progress and challenges. *Proceedings of the Combustion Institute*, 29 :1–28, 2002.

[8] N. Docquier et S. Candel. Combustion control and sensors : a review. *Progress in Energy and Combustion Science*, 28 :107–150, 2002.

[9] V. FAIVRE. *Etude expérimentale et numérique du contrôle actif des jets dans des chambres de combustion*. PhD thesis, Institut National Polytechnique de Toulouse, 2003.

[10] F. Bouziani, I. D. Landau, R. R. Bitmead, et A. Voda-Besançon. An analytically tractable model for combustion instability. *CDC-ECC*, (44), 2005.

[11] F. Bouziani, I. D. Landau, R. R. Bitmead, et A. Voda-Besançon. Analysis of two coupled Van der Pol equations as a model for combustion instability. Technical report, Laboratoire d'Automatique de Grenoble, ENSIEG BP 46, 38402 Saint-Martin d'Hères, France, Mar 2005.

[12] Wayne J. Dunstan. *System Identification of Nonlinear Resonant Systems*. PhD thesis, University of California, San Diego, 2003.

[13] W. J. Dunstan et R. R. Bitmead. Identification of resonant systems using periodic multiplicative reference systems. *13th IFAC Symposium on System Identification, Rotterdam NL*, 2003.

[14] A.A. Peracchio et W.M. Proscia. Nonlinear heat-release/acoustic model for thermoacoustic instability in lean premixed combustors. *ASME J. Eng. Gas Turbines Power*, 121(3):415-421, 1999.

[15] N.N. Bogoliubov et Y.A. Mitropolski. *Asymptotic Methods in the Theory of Nonlinear Oscillations*. Hindustan Publishing Corp, Delhi, and Gordon and Breach, New York, 1961.

[16] C. Hayashi. *Nonlinear Oscillations in Physical Systems*. McGraw-Hill Book Co, New York, 1964. (reprinted by Princeton University Press, 1985).

[17] P. S. Landa. *Nonlinear Oscillations and Waves in Dynamical Systems*. Kluwer, 1996.

[18] P. S. Landa. *Regular and Chaotic Oscillations*. Springer, New York, 2000.

[19] I. D. Landau et R. R. Bitmead. On the method of Krylov and Bogoliubov for the analysis of nonlinear oscillations. Technical report, Mechanical and Aerospace Engineering Department, University of California, San Diego, 9500 Gilman Drive, La Jolla CA 92093-0411, USA, Jan 2004.

[20] F. Bouziani, I. D. Landau, et A. Voda-Besançon. First and second-order K-B approximations for the analysis of nonlinear oscillations in autonomous systems. Technical report, Laboratoire d'Automatique de Grenoble, ENSIEG BP 46, 38402 Saint-Martin d'Hères, France, Jun 2004.