# Modélisation de l'enlèvement de matériau abradable à l'aide d'équations différentielles à retard

Nicolas Salvat, Alain Batailly, Mathias Legrand

► To cite this version:

Nicolas Salvat, Alain Batailly, Mathias Legrand. Modélisation de l'enlèvement de matériau abradable à l'aide d'équations différentielles à retard. 11e Colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01104675

HAL Id: hal-01104675

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

Submitted on 10 Nov 2016

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
# Modélisation de l'enlèvement de matériau abradable à l'aide d'équations différentielles à retard

Nicolas SALVAT, Alain BATAILLY, Mathias LEGRAND

Laboratoire de dynamique des structures et vibrations, département de génie mécanique, université McGill, Montréal, Québec

**Résumé** — L'amélioration de l'efficacité énergétique des turbomachines aéronautiques s'organise selon des axes de recherche distincts, l'un d'entre eux concerne la réduction des jeux entre parties tournantes et le carter environnant afin d'augmenter les taux de compression. La réduction des jeux a notamment pour effet de favoriser l'apparition de contacts structuraux qui sont dorénavant considérés comme des événements normaux, bien qu'indésirables, pouvant se produire au cours du fonctionnement de la turbomachine. Longtemps vus comme une solution robuste pour empêcher l'apparition de hauts niveaux de vibration lors de contacts structuraux, les matériaux abradables n'évitent toutefois pas certains scenarii d'interaction préjudiciables pour l'intégrité du moteur. Un examen qualitatif de ces phénomènes d'interaction est réalisé par l'intermédiaire d'une analyse de stabilité d'un système décrit à l'aide d'équations différentielles à retard. Les bases théoriques sont brièvement exposées puis les résultats sont comparés à ceux obtenus avec des approches plus classiques telles que des simulations temporelles. Les résultats présentés mettent en évidence le caractère prédictif de l'analyse de stabilité.

**Mots clés** — interaction rotor/stator, matériau abradable, équations différentielles à retard, analyse de stabilité

## 1 Introduction

Les équations différentielles à retard (EDR) se distinguent des équations différentielles ordinaires en ce sens que l'évolution du système à un instant $t$ dépend de l'histoire du système à $t - \tau$. Les termes de retard $\tau$ sont généralement des constantes positives ou des variables dépendant continûment du temps ou de l'état du système [1]. Dans le cadre de cette étude, seul un retard constant en déplacement est considéré :

$$ \frac{dz}{dt}(t) = \dot{z}(t) = g(t, z(t), z(t - \tau)) \quad (1) $$

Ce type d'équations nécessite des fonctions définies sur l'intervalle $[-\tau; 0]$ comme conditions initiales. Elles sont couramment utilisées en épidémiologie [2] ou dans le domaine de l'usinage [3] où elles permettent la modélisation des procédés de coupe afin d'optimiser les taux d'enlèvement de matière sans atteindre des zones d'instabilité vibratoire. Dans le cadre de notre étude, les EDR sont envisagées pour modéliser l'enlèvement de matériau abradable lors de contacts aube/carter dans des turbomachines aéronautiques. Ces revêtements abradables sont usuellement disposés sur le carter afin de minimiser le jeu de fonctionnement avec les sommets d'aubes tout en réduisant les efforts de contact en cas d'interaction. Ces matériaux doivent satisfaire un compromis mécanique : être suffisamment souples pour ne pas endommager les aubes lors de contacts mais également suffisamment résistants à leur environnement hostile (haute température, corrosion...). Longtemps entrevus comme une solution robuste pour empêcher l'apparition de hauts niveaux de vibrations, les matériaux abradables n'éliminent toutefois pas certains scénarii indésirables. Des travaux de natures expérimentale [4] et numérique [5] ont notamment mis en évidence des phénomènes d'interaction de type aube/abradable pouvant conduire à une rupture d'aube. Les premières études numériques envisagées ont favorisé des stratégies d'intégration temporelle explicite [6]. Toutefois les temps de calcul élevés motivent la mise en œuvre d'analyses purement qualitatives pour espérer une compréhension approfondie de ces phénomènes d'interaction.
Sous certaines conditions détaillées dans la section 3 une analogie entre ce phénomène d'enlèvement de matière et le fraisage est envisageable comme succinctement illustré sur la figure 1. Il s'agit notamment de supposer qu'une seule aube entre en contact avec l'abradable et qu'aucune séparation n'est possible. Le cadre mathématique des EDR utilisé dans le contexte du fraisage et du tournage est alors disponible et facilement exploitable pour explorer la dynamique de l'interaction envisagée. L'axe de rotation de l'arbre est parfaitement rigide et le contact est la conséquence de la vibration des aubes ou de la déformation du carter uniquement et l'équation du mouvement s'écrit :

$$M\ddot{u}(t) + C\dot{u}(t) + Ku(t) = f(t, u(t), u(t-\tau)) \quad (2)$$

où M, C et K sont respectivement les matrices masse, amortissement et raideur de l'aube étudiée, et f est le vecteur des efforts extérieurs dépendant de l'état du système aux temps t et $t - \tau$. Ce terme de retard provient de la rotation de l'aube et du fait que sa pénétration dans l'abradable à un temps t dépend du même événement au tour précédent.

<figure><img src="image_2.png" /><figcaption>Figure 1 - Analogie entre l'enlèvement d'abradable et le fraisage</figcaption></figure>

# 2 Équations différentielles à retard

## 2.1 Contexte théorique

Les EDR sont utiles pour la modélisation de systèmes comprenant des retards tels que l'équation (1). Dans le cas particulier des systèmes mécaniques, ce couplage est souvent pris en compte dans le terme d'effort f de l'équation (2) de la manière suivante :

$$f(t) = \Gamma(t)(u(t-\tau) - u(t)) \quad (3)$$

où $\Gamma(t)$ est la matrice des efforts extérieurs. L'introduction du vecteur d'état de dimension $2n$ $z^T(t) = (u^T(t) \ \dot{u}^T(t))$ permet de se ramener à une EDR du premier ordre :

$$\begin{aligned} \dot{z}(t) &= A(t)z(t) + B(t)z(t-\tau) \\ z(t) &= \varphi(t), \quad t \in [-\tau; 0] \end{aligned} \quad (4)$$

avec :

$$A(t) = \begin{bmatrix} 0_n & I_n \\ -M^{-1}(K + \Gamma(t)) & -M^{-1}C \end{bmatrix} \quad \text{et} \quad B(t) = \begin{bmatrix} 0_n & 0_n \\ M^{-1}\Gamma(t) & 0_n \end{bmatrix} \quad (5)$$

où $A(t)$ and $B(t)$ sont des matrices $2n \times 2n$ dépendantes des paramètres du système et des efforts extérieurs, et $\varphi$ est la fonction des conditions initiales définie sur l'intervalle $[-\tau; 0]$. $I_n$ et $0_n$ sont respectivement la matrice identité et la matrice nulle de dimension $n \times n$.

Deux cas sont d'un intérêt particulier : (1) les EDR à coefficients constants où $A(t) = A$ et $B(t) = B$, et (2) les EDR à coefficients périodiques de période $T$ où $A(t + T) = A(t)$ et $B(t + T) = B(t)^1$. Lorsque l'une de ces conditions est satisfaite, le système d'équations (4) et (5) se prête à réaliser des études de stabilité de la solution d'équilibre $z = 0$.

<footnote n="1">1. Dans ce papier les deux cas sont considérés et les développements théoriques se limitent au cas où $T = \tau$.</footnote>
## 2.2 Analyse de stabilité

Un premier type d'analyse pour des matrices A et B constantes est basé sur la résolution d'un problème aux valeurs propres non linéaire, en considérant une solution exponentielle — du type $z(t) = Z \exp(\lambda t)$ — afin d'obtenir l'équation caractéristique :

$$ (I_{2n}\lambda - A - B\exp(-\lambda\tau))Z = 0 \quad (6) $$

associée au système (4). Le retard introduit un terme exponentiel dans cette équation rendant le problème de dimension infinie [7]. Les méthodes de résolution sont généralement basées sur des approximations rationnelles ou polynomiales de cette équation, telles que la discrétisation de l'opérateur solution (SOD) et la discrétisation du générateur infinitésimal (IGD). L'ensemble infini dénombrable des $\lambda$ solutions de cette équation déterminent la stabilité de la solution d'équilibre $z(t) = 0$.

Un deuxième type d'analyse, qui s'étend aux matrices périodiques, suppose une discrétisation fine du temps sur l'intervalle $[-\tau, \tau]$ et repose sur la construction de la matrice de monodromie $\Phi$, reliant l'état du système au temps $t_i$ à celui de l'instant $t_i - \tau$, i.e. $z(t_i) = \Phi z(t_i - \tau)$. Cette matrice est une approximation finie de l'opérateur de monodromie qui agit sur l'espace des fonctions continues et qui est donc théoriquement de dimension infinie. Différentes méthodes existent au sein de cette famille, telles que les techniques de collocation, la méthode de semi-discrétisation (SDM), utilisée dans ce papier et détaillée dans l'annexe A, ou encore des analyses éléments finis en temps. Selon la théorie de Floquet, la stabilité de la solution $z(t) = 0$ est alors déterminée en calculant les valeurs propres de $\Phi$ [1].

Ces deux types d'analyse permettent la construction de diagrammes de stabilité dépendant de paramètres d'intérêt et comprenant des zones de fonctionnement connues sous le nom de lobes de stabilité, représentant la limite entre les zones stables et instables de la solution d'équilibre. Dans notre étude, ces deux paramètres sont la vitesse de rotation et le coefficient d'effort contenu dans $\Gamma(t)$. Dans ce type de stratégie, aucune intégration temporelle avec condition initiale n'est utile.

# 3 Modélisation

## 3.1 Aube

L'aube étudiée appartient au dernier étage du compresseur basse pression d'un moteur d'avion. Le maillage éléments finis contient 11 339 éléments (tétraèdres quadratiques) totalisant 22 898 nœuds.

Le contact entre le sommet d'aube et le carter est géré au bord d'attaque ou au bord de fuite uniquement. La très grande taille du modèle éléments finis ne permet pas d'envisager des temps de calculs raisonnables, il est donc réduit par une méthode de synthèse modale inspirée de la méthode de Craig-Bampton et permettant de prendre en compte les effets centrifuges [8]. L'équation du mouvement est projetée dans l'espace réduit :

$$ u = \begin{pmatrix} u_b \\ u_i \end{pmatrix} = \begin{bmatrix} I & 0 \\ \Psi_s(0) & \Psi_c \end{bmatrix} \begin{pmatrix} u_b \\ q \end{pmatrix} \quad (7) $$

Le modèle réduit est construit de façon à traiter le contact directement dans l'espace réduit sur les nœuds frontières $u_b$. Les déplacements intérieurs $u_i$ sont réduits à des participations modales $q$ de modes encastrés $\Psi_c$ calculés à différentes vitesses de rotation [9].

## 3.2 Loi de coupe en fraisage

Lors du fraisage, sous réserve que les vibrations entre l'outil et la pièce usinée sont négligeables, l'évolution de l'épaisseur du copeau $h$ est périodique et dépend de l'angle de coupe $\phi$ ainsi que de la profondeur de passe par dent $f_t$. L'effort de coupe $F$ est usuellement considéré proportionnel à l'épaisseur du copeau et est donc également périodique en temps :

$$ F(t) = K_e h(t) \quad \text{avec} \quad h(t) = f_t \sin(\phi(t)) \quad (8) $$

où $K_e$ est le coefficient d'effort spécifique et dépend du matériau de la pièce usinée et de la profondeur de passe [10].
<figure><img src="image_3.png" /><figcaption>Figure 2 - Analogie entre le fraisage et l'enlèvement d'abradable dans une turbomachine</figcaption></figure>

Lorsque du broutement se produit (du fait d'un couplage modal ou d'un état de surface ondulé) l'épaisseur du copeau varie en fonction de la position de l'outil par rapport à la pièce. La variation de l'effort de coupe qui en découle entraîne des vibrations auto induites [10]. La relation entre l'épaisseur du copeau $h$ et la position radiale de l'outil $u_r$ s'écrit alors :

$$h(t) = f_t \sin(\phi(t)) + u_r(t - \tau) - u_r(t) \quad (9)$$

ce qui est illustré par la figure 2(a).

La prise en compte des vibrations autour de la trajectoire nominale de la dent de l'outil se traduit par un système d'EDR couplées. Usuellement, l'espace d'état utilisé pour le fraisage ressemble à celui de l'équation (4) et associe un modèle à un ou deux degrés de liberté (ddl) de l'outil et/ou de la pièce usinée. L'étude de stabilité des solutions à ces équations est faite pour estimer le taux d'enlèvement de matière maximum possible [11].

## 3.3 EDR et enlèvement de matériau abradable

Dans cette section, seuls les déplacements radiaux des nœuds frontières sont pris en compte dans le vecteur $u_b$.

De même que pour le fraisage, un effort de coupe proportionnel à la pénétration de l'aube dans l'abradable est supposé pour représenter l'effort de contact en sommet d'aube. Cette pénétration est déterminée par l'évolution du déplacement du sommet d'aube entre deux tours consécutifs, comme le montre la figure 2(b), en supposant qu'il n'y ait jamais de séparation entre l'aube et le matériau abradable (hypothèse de contact permanent). L'effort de contact radial s'écrit alors :

$$f(h(t)) = K_e h(t) = K_e [u_b(t - \tau) - u_b(t)] \quad (10)$$

où l'épaisseur de matériau abradable enlevé est conservée en mémoire dans le vecteur $h$ qui dépend du déplacement des nœuds frontières $u_b$ d'un tour à l'autre et de la configuration de contact.

Les équations du mouvement sont similaires aux équations (2) et (3) où $\Gamma(t) = \Gamma$ est constante (du fait de l'hypothèse de contact permanent) et contient le coefficient de coupe $K_e$ appliqué sur les nœuds frontières sélectionnés. La gestion du contact est donc réalisée via une loi linéaire de type pénalité, sensée représenter de façon très simplifiée les efforts de contact calculés dans [5] à l'aide des multiplicateurs de Lagrange.

Une utilisation judicieuse des coordonnées modales $u = Vx$ permet de projeter les équations (2) et (3) :

$$\ddot{x}(t) + 2\xi\Lambda\dot{x}(t) + \Lambda^2 x(t) = V^T \Gamma V [x(t - \tau) - x(t)] \quad (11)$$

Ce système linéaire est remanié dans l'espace d'état pour obtenir une équation similaire à l'équation (4) où A et B sont des matrices constantes :

$$A = \begin{bmatrix} 0_n & I_n \\ -\Lambda^2 - V^T \Gamma V & -2\xi\Lambda \end{bmatrix} \quad \text{et} \quad B = \begin{bmatrix} 0_n & 0_n \\ V^T \Gamma V & 0_n \end{bmatrix} \quad (12)$$
## 3.4 Prise en compte du décollement

On présente brièvement une extension de la modélisation de l'enlèvement d'abradable par EDR avec prise en compte du décollement entre l'aube et le revêtement. L'hypothèse de contact permanent présentée dans la section précédente est une des limitations inhérente à l'utilisation des EDR linéaires mais est cependant obligatoire pour mener une étude de stabilité.

En effet, si on suppose que l'aube peut *se séparer* du revêtement abradable, l'éventuelle pénétration de l'aube dans l'abradable à l'instant $t$ ne peut plus être déduite de la différence $\mathbf{x}(t - \tau) - \mathbf{x}(t)$. Une variable auxiliaire $\chi(t)$ est introduite afin d'enregistrer le profil de la couche d'abradable et de tester à chaque instant si l'aube est en contact ou pas. Ainsi, il convient de distinguer deux états :

- *l'interaction* pendant laquelle l'aube pénètre dans le revêtement abradable, l'effort de coupe devenant fonction de $\chi(t - \tau)$ et conduisant à l'équation :

$$ \ddot{x}(t) + 2\xi\Lambda\dot{x}(t) + \Lambda^2 x(t) = V^T\Gamma V[\chi(t - \tau) - x(t)] \quad (13) $$

- *le décollement* au cours duquel le jeu aube/revêtement abradable est strictement positif et l'EDR (13) devient une équation différentielle classique sans terme de retard :

$$ \ddot{x}(t) + 2\xi\Lambda\dot{x}(t) + \Lambda^2 x(t) = 0 \quad (14) $$

Le problème devient non-linéaire et dépendant des conditions initiales. Toute analyse de stabilité comme établie précédemment devient caduque et les seules prédictions envisageables sont des intégrations temporelles au cours desquelles il est nécessaire d'alterner entre solveurs pour EDR et solveurs pour équations différentielles classiques. Ce cadre est un domaine à explorer mais dans ce travail, toutes les simulations supposent un contact permanent entre l'aube et le revêtement abradable, et emploient la méthode de semi-discrétisation présentée dans l'annexe A.

# 4 Résultats préliminaires

L'étude de stabilité par rapport au paramètre $K_e$ introduit dans l'équation (10), est conduite sur l'intervalle $K_e \in [0; 2, 5] \times 10^4$ N/m de telle sorte que les efforts de contact sont cohérents avec les valeurs rapportées dans [5]. Toutes les fréquences données dans cette section sont normées par rapport à la première fréquence propre de l'aube. La plage de vitesse de rotation considérée est $\Omega \in [0; 0, 45]$.

Les résultats sont présentés sous forme de lobes de stabilité, séparant les zones de fonctionnement où la solution $z = 0$ est stable, des zones où elle est instable. Pour rappel, la solution $z = 0$ correspond à une configuration dans laquelle l'aube tourne dans le carter avec un jeu aube/abradable nul sans vibrer.

## 4.1 Convergence en temps et en espace

Les lobes illustrés dans les figures 3 et 4 représentent la stabilité où les zones grises correspondent à des combinaisons $\Omega$; $K_e$ conduisant à des comportements divergents. Afin de s'assurer de la validité de ces lobes, deux types de convergence doivent être vérifiés. D'une part la convergence en temps, liée à la discrétisation de la période de révolution de l'aube, et d'autre part, la convergence en espace, associée au nombre de modes retenus. Dans un premier temps, les effets centrifuges sont négligés et le contact est uniquement géré au bord d'attaque. Le retard $\tau = 1/\Omega$ dépend de la vitesse de rotation $\Omega$ et est discrétisé en $N$ pas de temps $\Delta t : \tau = N\Delta t$. La convergence des simulations est vérifiée sur toute la plage de vitesse en augmentant $N$.

Les lobes sur la figure 3 mettent en évidence que pour des vitesses de rotation élevées ($\Omega > 0,25$), $N = 40$ est suffisant pour obtenir des résultats convergés. En revanche, pour de plus faibles vitesses, la convergence est atteinte uniquement pour $N \ge 60$. Par conséquent, la valeur $N = 60$ est utilisée dans toutes les simulations présentées. De manière similaire, la convergence en espace est vérifiée suivant le nombre de modes retenus dans l'équation (11). Les deux premiers modes de l'aube (flexion 1F et torsion 1T) sont suffisants.
<figure><img src="image_4.png" /><figcaption>**Figure 3** – Convergence en temps : N = 40 (—); N = 60 (---) et N = 80 (—-)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4** – Lobes d'instabilité avec (---) et sans (—) effets centrifuges</figcaption></figure>

## 4.2 Effets centrifuges

Les fréquences propres de l'aube sont modifiées du fait du raidissement centrifuge ce qui a un effet sur les lobes d'instabilité dessinés sur la figure 4. Il est remarquable que la prise en compte des effets centrifuges augmente le domaine de stabilité à haute vitesse. Les effets centrifuges sont pris en compte dans toutes les simulations suivantes.

<figure><img src="image_6.png" /><figcaption>**Figure 5** – Lobes d'instabilité : flexion (---) et torsion (—)</figcaption></figure>

# 5 Résultats

Dans cette section, l'étude de stabilité est menée à bien en ne conservant qu'une seule contribution modale. Deux ensembles de lobes sont superposés sur la figure 5, où il est visible que les lobes d'instabilité associés à la flexion sont plus larges et plus critiques (car se produisant pour de plus faibles valeurs de $K_e$ et pour des vitesses de rotation plus faibles) que ceux associés à la torsion.

La figure 6 contient le spectre du déplacement radial du bord d'attaque des simulations temporelles effectuées pour chaque vitesse de rotation en suivant la stratégie présentée dans [9]. Ces simulations temporelles sont réalisées en prenant en compte le contact uniquement au bord d'attaque, le reste des nœuds en sommet d'aube pouvant vibrer librement. Deux séries de pics sont visibles autour de chacune des deux premières fréquences propres de l'aube.

Le long de la ligne en pointillés bleus sur la figure 6, la position des pics est en accord avec les lobes d'instabilité associés au mode de flexion. Cette observation est également valable pour le premier mode de torsion et met en valeur la cohérence entre les résultats obtenus par simulations temporelles et par analyses de stabilité.

Les résultats obtenus avec les EDR révèlent que l'aube est plus sensible à l'instabilité en flexion qu'en torsion. Ces résultats sont confirmés par les plus grands pics observés autour de la première fréquence propre sur les spectres des simulations temporelles et sont en accord avec des travaux précédents [5] où le rôle essentiel du premier mode de flexion est mis en avant.

Les pics prédits par les simulations temporelles sont en accord avec les lobes d'instabilité des EDR malgré les hypothèses contraignantes, notamment le contact permanent aube/abradable. Le fait que deux modes suffisent à prédire les zones d'instabilité est également en accord avec les résultats temporels puisque les spectres obtenus ne présentent pas de pics significatifs pour de hautes fréquences $f > 5$ non représentées sur la figure 6.
<figure><img src="image_7.png" /><figcaption>**Figure 6** – Spectre des simulations temporelles : premier mode de torsion (—) et premier mode de flexion (---); les lobes d'instabilité de la figure 5 sont superposés.</figcaption></figure>

Enfin, il semble intéressant de remarquer que les pics d'amplitude détectés pour de faibles vitesses de rotation, dont l'existence était supposée liée au contact unilatéral [5, 6], sont prédits par une EDR sans terme non linéaire. Ces résultats indiquent que l'amplitude de ces pics est probablement autant dépendante de la nature régénérative du phénomène d'enlèvement de matière qu'à la non linéarité induite par le contact unilatéral.

# 6 Conclusion

L'enlèvement de revêtement abreadable par les aubes dans un moteur d'avion représente un problème complexe, jusqu'à présent uniquement traité par des techniques d'intégration temporelle. Dans le cadre de ce travail, il est mis en évidence que l'utilisation d'Equations Différentielles à Retard (EDR) linéaires permet de mener des analyses de stabilité directes issues d'un problème aux valeurs propres. Les résultats confirment la pertinence de cette stratégie puisqu'ils confirment l'approche temporelle existante.

Les premiers résultats illustrent la convergence de la méthode et la précision des lobes d'instabilité et ce, en ne conservant qu'un faible nombre de modes d'aube. La comparaison avec les simulations temporelles souligne la cohérence des deux approches. Les EDR permettent de prédire précisément les pics de réponse de l'aube. Qui plus est, l'analyse de stabilité réduit significativement les temps de calculs. Elle pourrait être incorporée au processus d'optimisation du profil d'une aube dès la phase de conception pour maximiser le domaine de stabilité au sein d'une plage de vitesse donnée. Les concepteurs pourraient bénéficier de cette procédure pour éliminer les profils d'aube peu robustes aux interactions aube/abreadable.

Parmi les futurs développements envisageables, des travaux sont en cours pour évaluer (1) l'importance du frottement dans le phénomène d'interaction, (2) les effets de la séparation aube/revêtement abreadable et (3) le couplage avec les déplacements tangentiels en sommet d'aube. Enfin, l'implémentation d'un carter flexible nécessitant des EDR couplées est également en cours de développement.

# A Méthode de Semi-Discrétisation

La méthode de semi-discrétisation (SDM) a récemment été développée pour l'analyse de stabilité des EDR linéaires périodiques en temps. Elle repose sur le calcul des valeurs propres de l'opérateur de monodromie. La théorie de Floquet est ensuite utilisée pour l'analyse de stabilité. La SDM a notamment été utilisée pour l'étude du tournage [12].

Pour une EDR périodique en temps de dimension $2n$ telle que l'équation (4), la SDM repose sur la discrétisation du terme de retard $z(t - \tau)$ en une fonction constante par morceaux sur une période $\tau = N\Delta t$ alors que tous les autres termes dépendant du temps restent inchangés. Pour $\Delta t$ suffisamment petit (ou $N$ suffisamment grand), on a $\Delta t = [t_i ; t_{i+1}]$, avec $t_i \in [0 ; \tau]$ et $i = 0, 1, \dots (N - 1)$. Ainsi,
∀t ∈ [tᵢ; tᵢ₊₁], l'équation (4) est équivalente à :

$$ \dot{z}(t) = A(t)z(t) + B(t)z_{\tau,i} \quad (15) $$

où $z_{\tau,i}$ est une combinaison linéaire pondérée des vecteurs retardés $z(t_{i-N})$ et $z(t_{i-N+1})$ :

$$ z_{\tau,i} = \frac{1}{2}z_{i-N+1} + \frac{1}{2}z_{i-N} \approx z(t-\tau) \quad (16) $$

Le terme retardé est ainsi supposé constant sur l'intervalle $[t_i ; t_{i+1}]$ permettant ainsi une approximation de l'équation (4) par une équation différentielle classique de dimension finie pour laquelle une solution analytique existe :

$$ z(t) = \exp(A_i(t - t_i))(z_i + A_i^{-1}B_i z_{\tau,i}) - A_i^{-1}B_i z_{\tau,i} \quad (17) $$

où $z(t_i) = z_i$ est une condition initiale déterminée, et les coefficients linéaires discrets sont approchés par leur valeur moyenne sur l'intervalle choisi :

$$ A_i = \frac{1}{\Delta t} \int_{t_i}^{t_{i+1}} A(t)dt \quad \text{et} \quad B_i = \frac{1}{\Delta t} \int_{t_i}^{t_{i+1}} B(t)dt \quad (18) $$

La combinaison des équations (16) et (17) pour $t = t_{i+1}$ entraîne :

$$ z_{i+1} = P_i z_i + \frac{1}{2} R_i z_{i-N+1} + \frac{1}{2} R_i z_{i-N} \quad (19) $$

avec $P_i = \exp(A_i \Delta t)$ et $R_i = (\exp(A_i \Delta t) - I_{2n})A_i^{-1}B_i$. Une carte discrète peut ensuite être créée pour définir un vecteur d'état colonne $y_i$ de dimension $n(N+1)$ :

$$ y_i = [u_i \quad \dot{u}_i \quad u_{i-1} \quad u_{i-2} \quad \dots \quad u_{i-N}]^T \quad \text{et} \quad y_{i+1} = D_i y_i \quad (20) $$

avec :

$$ D_i = \begin{bmatrix} P_{i,1:n} & P_{i,n+1:2n} & 0_n & \cdots & 0_n & \frac{1}{2}R_{i,1:n} & \frac{1}{2}R_{i,n+1:n} \\ I_n & 0_n & 0_n & \cdots & 0_n & 0_n & 0_n \\ 0_n & 0_n & I_n & & 0_n & 0_n & 0_n \\ \vdots & & \ddots & & \ddots & & \vdots \\ 0_n & 0_n & & I_n & 0_n & 0_n & 0_n \\ 0_n & 0_n & \cdots & 0_n & I_n & 0_n & 0_n \end{bmatrix} \quad (21) $$

où $P_{i,1:n}$ représente les colonnes 1 à n de la matrice $P_i$. Comme le retard s'applique uniquement aux déplacements, seules les n premières colonnes de la matrice $R_i$ sont prises en compte et les termes retardés en vitesse ne sont pas inclus. Une approximation de l'opérateur de monodromie $\Phi$ peut ensuite être construite :

$$ y_N = \Phi y_0 \quad \text{avec} \quad \Phi = D_{N-1}D_{N-2}\cdots D_1D_0 \quad (22) $$

Le nombre de multiplications requises pour construire $\Phi$ dépend uniquement du paramètre N alors que la taille du système dépend à la fois de la discrétisation spatiale et temporelle. Finalement, l'étude de stabilité se concentre sur l'étude du module des valeurs propres de $\Phi$.

## Références

[1] B. Balachandran, T. Kalmar-Nagy, and D. Gilsinn. *Delay differential equations: Recent advances and new directions*. Springer Verlag, 2009.

[2] K. Cooke, P. Van den Driessche, and X. Zou. Interaction of maturation delay and nonlinear birth in population and epidemic models. *Journal of Mathematical biology*, 39:332–352, 1999.

[3] G. Quintana and J. Ciurana. Chatter in machining processes: A review. *International Journal of Machine Tools and Manufacture*, 51:363–376, 2011.
[4] A. Millecamps, J.-F. Brunel, P. Dufrénoy, F. Garcin, and M. Nucci. Influence of thermal effects during blade-casing contact experiments. In *ASME Conference Proceedings: 22nd Biennial Conference on Mechanical Vibration and Noise*, number 48982, pages 855–862, San Diego, USA, 2009.

[5] A. Batailly, M. Legrand, A. Millecamps, and F. Garcin. Numerical-experimental comparison in the simulation of rotor/stator interaction through blade-tip/abradable coating contact. *Journal of Engineering for Gas Turbines and Power*, 134(8), 2012.

[6] M. Legrand, A. Batailly, B. Magnain, P. Cartraud, and C. Pierre. Full three dimensional investigation of structural contact interactions in turbomachines. *Journal of Sound and Vibration*, 331:2578–2601, 2012.

[7] E. Jarlebring. *The spectrum of delay-differential equations: numerical methods, stability and perturbation*. PhD thesis, Technische Universität Carolo-Wilhelmina zu Braunschweig, Vännäs, Schweden, 2008.

[8] A. Sternchüss and E. Balmès. On the reduction of quasi-cyclic disks with variable rotation speeds. *Proceedings of the International Conference on Advanced Acoustics and Vibration Engineering (ISMA)*, pages 3925–3939, 2006.

[9] M. Legrand, A. Batailly, and C. Pierre. Numerical investigation of abradable coating removal through plastic constitutive law in aircraft engine. *Journal of Computational and Nonlinear Dynamics*, 7, 2012.

[10] J. Tlusty. *Manufacturing processes and equipment*. Prentice Hall, 2000.

[11] G. Stépán, R. Szalai, B. Mann, P. Bayly, T. Insperger, J. Gradisek, and E. Govekar. Nonlinear dynamics of high-speed milling—analyses, numerics, and experiments. *Journal of Vibration and Acoustics*, 127(2):197–203, 2005.

[12] T. Insperger and G. Stépán. *Semi-discretization for time-delay systems: stability and engineering applications*. Springer Science, 2011.