FR 3089266 A1 20200605 FR 1872222 A 20181203 Description Titre de l'invention : PRÉDICTION D'INSTABILITÉ DE CRISSEMENT POUR UN FREIN DE VÉHICULE Domaine technique [0001] La présente invention concerne une méthode et un dispositif de prédiction d'instabilité de crissement pour un système de frein, présentant au moins une surface de friction, par exemple un frein de véhicule à disque, à tambour, etc. [0002] Il est difficile de prédire l'apparition d'un bruit de crissement suffisamment tôt pour réduire les coûts de développement dans le processus de fabrication d'un système de freinage. Par conséquent le crissement est toujours responsable d'un nombre important de nuisances dans le domaine de l'automobile et le problème de la prédiction du bruit de crissement et de la compréhension des mécanismes (instabilités) à son origine sont actuellement un enjeu pour les industriels de l'automobile. Technique antérieure [0003] La littérature recense de nombreuses méthodes pour l'étude des systèmes dy- namiques non-linéaires. Le choix de la méthode se fait selon certains critères tels que le type de non-linéarités impliquées ou encore le type de solution recherchée. En effet, ce ne sont pas les mêmes outils et les mêmes formalismes qui sont utilisés pour traiter des non-linéarités régulières ou non, fortes ou faibles, avec ou sans hystérésis et des solutions périodique, quasi-périodique, transitoire ou chaotique. [0004] Dans le cadre du crissement, il est nécessaire d'utiliser des méthodes permettant de traiter des systèmes autonomes non-linéaires dont les réponses sont souvent considérées comme étant périodiques ou quasi-périodiques. La taille des systèmes étudiés, les non-linéarités impliquées et les différents couplages rendent la résolution analytique des équations associées difficilement inenvisageable. [0005] Parmi les méthodes numériques en mesure d'tuudyser ces systèmes, certaines permettent cependant de statuer sur la stabilité de la position d'équilibre, comme par exemple, la méthode CEA (« Complex Eigenvalue Analysis ») qui met en oeuvre une analyse en valeur propre classique au voisinage de l'équilibre statique non-linéaire. Cette méthode est décrite de manière générale par exemple dans le document de MOIROT, Franck. intitulé « Etude de la stabilité d'un équilibre en présence de frottement de Coulomb. Application au crissement des freins à disque ». (1998, Thèse de doctorat, Palaiseau, Ecole polytechnique). [0006] La méthode CEA consiste à déterminer les valeurs propres du système autour (le sa position d'équilibre afin de statuer sur la stabilité de ce point d'équilibre. Les modes qui possèdent une valeur propre à partie réelle positive sont considérés comme 2 instables puisqu'une légère perturbation de l'équilibre entrainera un taux d'accroissement positif des oscillations associées à ces modes. [0007] Cette méthode CEA est relativement simple à mettre en oeuvre et ses résultats sont faciles à analyser puisque les modes instables identifiés sont potentiellement les modes instables présents dans la réponse dynamique du système. Cependant, la méthode CEA est souvent sur- ou sous-prédictive, c'est-à-dire que l'ensemble des modes instables identifiés par la méthode CEA ne sont pas nécessairement présents dans la réponse dynamique et que la réponse dynamique peut comprendre des modes instables (i.e. fréquences) non identifiées par la méthode CEA, correspondant à l'émergence de nouvelles instabilités. [0008] Il apparaît ainsi un besoin de trouver un critère et/ou une méthode permettant d'identifier dans la réponse dynamique du système de frein les modes instables qui sont responsables du crissement à partir des résultats de la méthode CEA. [0009] Par ailleurs, la méthode CEA permet uniquement d'identifier des modes poten- tiellement impliqués dans la réponse dynamique, mais elle ne donne aucune indication sur les niveaux vibratoires et les fréquences d'excitation du système. [0010] Ainsi l'obtention de la réponse dynamique et des niveaux vibratoires associés doit toutefois se faire avec d'autres méthodes numériques qui viennent compléter la méthode CEA et qui permettent d'approximer la réponse dynamique, en termes notamment de niveaux vibratoires et de fréquences d'excitation du système. [0011] La méthode MASA (Modal Amplitude Stability Analysis) a pour objectif de fournir une analyse complémentaire à la CEA. La méthode MASA est décrite de manière générale dans le document de NACIVET, S. et J.-J. SINOU. 2017, intitule «Modal amplitude stability analysis and its application to brake squeal», Applied Acoustics, vol. 116, p. 127-138. [0012] La MASA permet de déterminer les occurrences de crissement en estimant une ap- proximation de la solution dynamique initiale linéarisée autour d'un point d'équilibre pour une amplitude modale fixée ainsi que la contribution des différents modes potentiellement instables en évaluant la montée en divergence ou l'atténuation des modes potentiellement instables. L'objectif de la MASA est ainsi de fournir une analyse complémentaire à la CEA en identifiant les modes instables du système de frein. [0013] Un exemple d'application de cette méthode MASA aux instabilités de crissement d'un système de frein est décrit dans le document brevet publié sous le numéro FR3039645B1. [0014] Ce document brevet FR3039645B1 divulgue une méthode d'estimation des oc- currences de crissement dans des systèmes de frein qui permet d'estimer non seulement la stabilité du point d'équilibre non linéaire mais également l'apparition de fréquences instables causées par la génération du bruit de crissement. 3 100151 Cette méthode est basée sur un modèle numérique non linéaire du système de frein et permet d'obtenir une approximation de l'évolution clans le temps des parties réelles de valeurs propres complexes quand les amplitudes des vibrations non linéaires initiales augmentent. [0016] Cependant, dans la méthode décrite dans cc document brevet, l'analyse des modes instables est effectuée mode par mode, sur la base des modes identifiés par exemple une méthode telle que la méthode CEA. Une telle approche mono-mode permet d'identifier d'éventuels modes instables que la méthode CEA n'aurait pu identifier. Pour autant, elle ne permet pas de prendre en compte leurs participations. [0017] Cette approche mono-mode ne permet pas non plus de prendre en compte les couplages pouvant apparaitre entre les différentes instabilités et l'impact que ces couplages peuvent avoir sur les participations de chaque mode. [0018] Le brevet FR2988503B1 décrit un procédé de modélisation en vibration d'un système mécanique, tel un système de freinage, en un modèle réduit formé par des super-éléments. Résumé de l'invention [0019] Selon un premier aspect, l'invention concerne un procédé de prédiction d'instabilités de crissement d'une surface de friction d'un système de frein, le procédé étant mis en oeuvre par ordinateur et comprenant : a) une identification de N modes instables autour d'un point d'équilibre statique, les modes instables étant définis respectivement par N fréquences et N valeurs propres initiales de N modes du système de frein ; b) une modélisation du système de frein par un modèle linéaire dans lequel les efforts non-linéaires sont linéarisés vis-à-vis du déplacement et de la vitesse selon des degrés de liberté du système du frein ; c) une simulation d'une excitation du système de frein en imposant un déplacement selon lesdits degrés de liberté, le déplacement imposé étant sous forme de combinaison linéaire de N déformées modales, une déformée modale correspondant à un mode instable et étant définie par une amplitude modale de niveau variable et une fréquence parmi les N fréquences; d) un calcul d'une évolution d'au moins une valeur propre complexe d'au moins une matrice d'état du modèle linéaire lorsque le niveau d'une amplitude modale d'au moins une déformée modale varie; et c) une identification d'instabilités de crissement du système de frein sur la base de valeurs propres complexes calculées pour les niveaux d'amplitudes modales considérés. [0020] La nouvelle méthode proposée ici permet de suivre, au travers de l'évolution en fonction du temps des valeurs propres, l'évolution des modes instables identifiés initialement par exemple une méthode telle que la méthode CEA et de détecter si ces modes sont effectivement instables en régime dynamique. Cette méthode permet également d'identifier, au travers de l'évolution en fonction du temps des valeurs propres, d'éventuels modes instables que la méthode CEA n'aurait pu identifier. [0021] Selon au moins un mode de réalisation du procédé selon le premier aspect, une valeur courante d'une amplitude modale est définie en fonction du temps par récurrence par rapport à une valeur courante précédente et une amplitude modale initiale, et le procédé comprend un calcul de l'évolution en fonction du temps des valeurs propres complexes pour les niveaux d'amplitudes modales considérés. [0022] Selon au moins un mode de réalisation du procédé selon le premier aspect, comprenant le calcul de la partie réelle de chaque valeur propre complexe obtenue et de la fréquence correspondante. [0023] Selon au moins un mode de réalisation du procédé selon le premier aspect, dans lequel une valeur courante pk(LI) à l'instant t,, = t' + bt d'une amplitude modale d'un mode k est définie par rapport à la valeur précédente pk(t') à l'instant t' et à la valeur propre 21n k [0024] obtenue pour le mode k à l'itération précédente de la façon suivante [0025] [Math.101] [0026] Selon au moins un mode de réalisation du procédé selon le premier aspect, dans lequel le système de frein est décomposé en sous-systèmes dynamiques et une matrice d'état du modèle linéaire correspond à une instabilité et un sous-système dynamique du système de frein. [0027] Selon au moins un mode de réalisation du procédé selon le premier aspect, le procédé comprenant en outre : une interruption du calcul de l'évolution des valeurs propres du modèle linéaire sur détermination qu'une nouvelle valeur propre obtenue à l'étape d) ne fait pas partie des N valeurs propres initiales du système de frein et que la partie réelle de la nouvelle valeur propre est positive ; une modification de la modélisation des déplacements de sorte à introduire une fonction de déformation supplémentaire associée à la nouvelle valeur propre et à représenter les déplacements comme combinaison linéaire des N fonctions de déformation et de la fonction de déformation supplémentaire ; une reprise de l'étape d) de calcul de l'évolution des valeurs propres du modèle linéaire sur la base de la modélisation modifiée des déplacements. [0028] Selon un deuxième aspect, l'invention concerne un procédé d'optimisation d'un système de frein, comprenant : une obtention de valeurs de paramètres de conception structurelle du système de frein ; une génération d'un modèle linaire du système de frein à partir des valeurs obtenues ; une identification d'instabilités de crissement du système de frein en mettant en oeuvre un procédé de prédiction d'instabilités de crissement selon l'une quelconque des revendications précédentes sur la base du modèle linéaire généré ; une correction d'au moins une valeur d'au moins un paramètre de conception structurelle du système de frein sur la base des instabilités de crissement identifiées ; une correction du modèle linaire du système de frein à partir de ladite au moins une valeur corrigée ; une identification d'instabilités de crissement du système de frein en mettant en oeuvre un procédé de prédiction d'instabilités de crissement selon l'une quelconque des revendications précédentes sur la base du modèle linéaire corrigé. [0029] Selon un troisième aspect, l'invention concerne un dispositif comprenant des moyens de mise en oeuvre des étapes d'un procédé selon le premier ou deuxième aspect. [0030] Selon un quatrième aspect, l'invention concerne un programme d'ordinateur comportant des instructions de programme pour l'exécution des étapes d'un procédé selon le premier ou deuxième aspect. [0031] Selon un quatrième aspect, l'invention concerne un support d'enregistrement lisible par un ordinateur sur lequel est enregistré un programme d'ordinateur comprenant des instructions de programme pour l'exécution des étapes d'un procédé selon le premier ou deuxième aspect. Brève description des dessins [0032] D'autres avantages et exemples de réalisations sont décrits ci-après par référence aux figures dans lesquelles : [0033] [fig.1] représente un exemple dc système de frein. [0034] [fig.2] est un organigramme d'un procédé dc prédiction d'instabilité de crissement pour un système de frein selon la méthode MASA mono-mode ; [0035] [fig.3] est un exemple dc fonction utilisée dans un procédé dc prédiction d'instabilité dc crissement pour un système de frein selon un exemple de réalisation ; [0036] [fig.4] illustre des aspects mathématiques d'un procédé de prédiction d'instabilité de crissement pour un système de frein selon un exemple de réalisation ; [0037] [fig.5] est un organigramme d'un procédé de prédiction d'instabilité de crissement pour un système de frein selon un exemple de réalisation qui applique la méthode MASA multi-mode ; [0038] [fig.6A] 6 présente un exemple de courbes d'évolution en fonction de l'amplitude modale de parties réelles des valeurs propres correspondant à différentes modes instables d'un système de frein selon un exemple de réalisation ; [0039] [fig.6B] présente des exemples de courbes d'évolution en fonction de l'amplitude modale de fréquences correspondant à différentes modes instables d'un système de frein selon l'exemple de de la Fig. 6A; [0040] [fig.7A] présente un exemple de courbes d'évolution temporelle d'amplitude modale, de parties réelles des valeurs propres et fréquences correspondant à différentes modes instables d'un système de frein selon un exemple de réalisation ; [0041] [fig.7B] présente un exemple de courbes d'évolution temporelle de parties réelles de valeurs propres selon l'exemple de la Fig. 7A ; [0042] [fig.7C] présente un exemple de courbes d'évolution temporelle de parties réelles de fréquences selon l'exemple de la Fig. 7A ; [0043] [fig.8A] présente des courbes d'évolution temporelle de parties réelles des valeurs propres correspondant à un premier mode instable d'un système de frein selon un exemple de réalisation ; [0044] [fig.8B] présente des courbes d'évolution temporelle de parties réelles des valeurs de fréquences correspondant à un premier mode instable selon l'exemple de la Fig. 8A ; [0045] [fig.8C] présente des courbes d'évolution temporelle de parties réelles des valeurs propres correspondant à un deuxième mode instable d'un système de frein selon un exemple de réalisation ; [0046] [fig.8D] présente des courbes d'évolution temporelle de parties réelles des valeurs de fréquences correspondant à un deuxième mode instable selon l'exemple de la Fig. 8C . Description des modes de réalisation [0047] Dans ce document, un modèle en éléments finis d'un système de frein automobile in- dustriel est décrit. [0048] Un exemple de système de frein est présenté en FIG. 1. Ce système de frein est par exemple composé d'une roue 10, d'un disque 20, d'un moyeu 30, d'une articulation 40, d'un étrier 50, d'un piston 60, d'un support 70, d'une plaquette interne 80 et d'une 5 7 plaquette externe 90. Les tampons 80 et 90 appuient sur le disque 20, en fonction de l'intensité de la commande. [0049] Dans un premier temps, on rappelle ici le formalisme de la méthode MASA dans le cas de l'analyse mono-mode. La méthode MASA mono-mode repose sur une approximation du système non-linéaire par décomposition sur le premier harmonique et une nouvelle linéarisation des forces non linéaires. Les valeurs propres sont obtenues en effectuant une analyse de stabilité par rapport à l'amplitude modale. Une évolution des valeurs propres en fonction de l'amplitude modale est obtenue. [0050] L'équation dynamique du système de frein peut s'écrire sous la forme : [0051] [Math.1] ± + KX 1:,X1eip(X) + Fetp(X) + Xe [0052] dans laquelle X est. le vecteur des déplacements généralisés et où le point représente la dérivée par rapport au temps. M, C et K sont respectivement les matrices de masse, d'amortissement et de raideur structurelles du système. F e,e sont. les efforts extérieurs appliqués au système, à savoir la pression appliquée au piston et à l'étrier. Cela correspond à une action par un conducteur du véhicule sur la pédale de frein. F , F cl et F dp sont les contributions apportées par les interfaces des cieux sous-systèmes considérés où pour chaque interface le décollement est autorisé. Enfin, F id représente les efforts non-linéaires de contact et de frottement à l'interface de friction, par exemple l'interface entre plaquettes et disque dans le cas d'un frein à disque. [0053] En pratique, le système de frein est chargé de la manière suivante : une force nodale est distribuée uniformément sur la surface de piston et son opposée est distribuée sur l'étrier de la même manière. Comme expliqué précédemment dans l'article Fazio et al, Applied Acoustics, volume 91, 15 avril 2015, pages 12 à 24, les non-linéarités aux interfaces de friction comprennent à la fois les termes non-linéaires cubiques et la potentielle perte de contact entre les surfaces du disque et les surfaces de plaquette. Un champ de vitesse correspondant à la rotation du disque est imposé sur la surface du disque. De cette manière les forces de frictions sont déduites des forces de contact normales en utilisant la loi de Coulomb classique. Un état de glissement permanent est envisagé et on fait l'hypothèse d'un coefficient de friction constant. [0054] Analyse de stabilité [0055] L'analyse de stabilité peut être effectuée de la manière suivante. Une position d'équilibre statique U s est déterminée en résolvant l'équation - [0056] [Math.2] 1W5 + Fpip Ris) f&pall.$) Fcip ais) + [Us.] = Pen [0057] On définit une matrice K ', , T,, comme suit : 100581 Inah.31 8 KnIbtifs 100591 où les J k sont les matrices Jacobiennes des efforts F k au point d'équilibre statique U s . où k identifie l'une des interfaces parmi les interfaces désignées par p/p, e/p et c/p, avec p/p = piston/plaquette, e/p = étrier/plaquette et c/p = chape/plaquette. 100601 Ainsi, en considérant une légère perturbation AU autour du point d'équilibre statique U s , l'équation Math 1 devient : 100611 IMath.41 MAÙ ± ± (1C+ Kni.,LOAU ± AU) - I [0062] Le développement de Taylor au premier ordre des efforts non-linéaires autour de la position d'équilibre donne : [0063] [Math.5] 'fie (Us + .M1) (Us) + LIAI/ + [0064] et le problème aux valeurs propres associé est : [0065] [Math.6] M ± ± (I( Kni,Us j'il = [0066] Ces développements correspondent à la méthode CEA classique et permettent de dé- terminer les valeurs propres complexes ?,/ et leur mode propre associé 'loi La stabilité du système est déduite en fonction de la partie réelle des valeurs propres 2,/ . Nouvelle linéarisation du problème [0067] A partir de l'équation Math 5, on pose 52 tel que : [0068] [Math.7] Foi (Us + AU) = F t'Us) [0069] où n représente l'erreur commise sur le développement de Taylor à l'ordre 1. Ainsi, l'équation Math 4 peut se réécrire : 100701 IMath.81 MAU ± CAU K+Knins +-LOAU + fi = 0 [0071] Lorsque SL = 9, le problème aux valeurs propres associé à l'équation Math 8 correspond à l'équation Math 6 et permet d'assurer la continuité de la méthode. [0072] La solution recherchée du problème est une solution périodique, elle peut donc être approximée par une série de Fourier tronquée à l'ordre N h . En d'autres termes, la perturbation AU peut s'écrire : [0073] [Math.9] 9 AMI.) z.zz - Alla.{;1 tu-20 [0074] où AU k (I ) est la k'è" harmonique de AU(t) dans le domaine temporel et de coef- ficients de Fourier associés AUc et AUs. w correspond à la pulsation des cycles limites non-linéaires. [0075] De la même manière, les efforts non-linéaires !L peuvent être décomposés sous la forme d'une série de Fourier tronquée : [0076] [Math.10] N, i2 S2 t) k [0077] L'expression des efforts non-linéaires F ', en fonction du temps est connue de façon analytique. Ainsi, les efforts Sà(t) peuvent être exprimés dans le domaine temporel et l'utilisation d'un algorithme tel que l'algorithme AFT [Alternate Frequency Time], alternant les domaines fréquenticl et temporel, permet de déterminer les coefficients de Fourier e' et :4: [0078] L'algorithme AFT est décrit par exemple dans le document CAMERON, T. M. et GRIFFIN, J. H. An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems. Journal of applied mechanics, 1989, vol. 56, no 1, p. 149-154 [0079] La méthode part du postulat que le caractère stable ou instable d'un mode est lié à sa première harmonique et que donc l'étude de la seule stabilité d'un système peut se faire en se limitant à une unique harmonique. L'injection des équations Math 9 et Math 10 dans l'équation Math 8 et la troncature à la première harmonique donne le système dynamique suivant : [0080] [Math.11] ii ± CA ti ni) U j =O [0081] où : [0082] [Math.12] Vo- - - - sta: : - tu- 100831 et [0084] IMath.131 n sko o [0085] L'originalité de la méthode réside dans la linéarisation de n sous la forme : [0086] [Math.14] fli =KiAU C [0087] Cette linéarisation composée d'un terme de raideur équivalente et d'amortissement équivalent permet de considérer le déphasage entre les composantes de chargement de Q, et les déplacements AU,. De plus, les matrices KI et CI peuvent être déterminées analytiquement. La formulation analytique est donnée ultérieurement. Stabilité du nouveau système 100881 En injectant l'équation Math 14 dans l'équation Math Il, il vient : 100891 IMath.l5l MA fil (C (K + Knus +1,1 + KI ) A = 0 [0090] qui peut être réécrit en utilisant la variable d'état Y [0091] [Math.16] Atjt 100921 sous la forme : 100931 IMath.171 = [0094] avec la matrice dynamique A : [0095] [Math.l8] î (K+ KnUth "t /rat +-Ki) C) [0096] Conformément à la théorie de la stabilité des systèmes dynamiques, la stabilité du système selon l'équation Math 16 est déterminée en calculant les valeurs propres de la matrice A. Ainsi, la partie réelle associée à la valeur propre du mode instable j permet de statuer sur la divergence du mode j, un mode instable correspondant à une partie réelle positive. En outre, la partie imaginaire associée correspond à la pulsation du mode/ qui participe à la dynamique de la solution linéarisée autour du point d'équilibre pour une amplitude modale initiale fixée. Si d'autres valeurs propres ont une partie réelle qui devient positive, alors il y a un potentiel risque de divergence des modes correspondants Aussi il est pertinent d'utiliser dans ces conditions une MASA avec une approche multi-fréquentielle afin de prendre en considération ces modes dans la solution approchée et de statuer sur leur divergence. Cette méthode MASA multi- 11 fréquentielle (appelée aussi multi-mode) sera décrite plus en détail ci-dessous. Expression des déplacements 100971 Les déplacements AU correspondent aux déplacements des degrés de liberté du système de frein. 100981 Le principe de la méthode MASA consiste à étudier la stabilité du système pour des niveaux de vibration donnés. C'est-à-dire que les déplacements AU sont fixés (i.e. imposés par la simulation ) et que les efforts non linéaires Q, et donc la matrice d'état dynamique A, sont déterminés en fonction des déplacements AU imposés. La linéarisation et l'analyse dc stabilité faite sur la matrice d'état dynamique A sont valables uniquement pour les niveaux de vibration considérés. L'écriture des déplacements AU est expliquée ci-après. [0099] Supposons que le/ ème mode soit instable et dc valeur propre k j= a + alors l'écriture Math 12 est équivalente à l'écriture suivante : [0100] [Math.19] 101011 où 101021 IMath.201 frs.lf,- A 1.1- = [0103] où X, représente la déformée modale non normalisée. [0104] On notera au passage que l'écriture de la solution à la première harmonique introduit naturellement la condition de stabilité donnée par [0105] [Math.21] [0106] Si les déplacements AU sont non nuls et connus, alors ki est déduit de la matrice A en résolvant : [0107] [Math.22] [0108] avec : [0109] ]Math.23] [0110] Où [0111] 12 [Math.24] [0112] La reponse périodique d'un tel système est principalement pilotée par le mode instable identifié linéarisé autour de la position d'équilibre, c'est-à-dire que l'on peut approximer AU(t) par [0113] IMath.251 coiffe/0)ft ±iJe tait) [0114] Où c est un coefficient réel et 141 la déformée (appelée aussi déformée modale) du mode instable j. L'idée ici consiste à considérer différents niveaux d'amplitude de ce mode et pour chacun de ces niveaux d'amplitude à statuer sur la stabilité de la solution linéarisée pour le niveau d'amplitude considéré. [0115] Ainsi, si un mode instable de mode propre -140 a été identifié au cours de la CEA avec une valeur propre associée 10116] IMath.261 [01 17] où a; est la partie réelle supposée positive ici et ce sa pulsation, alors on considère différents niveaux d'amplitude modale J. - , où m fait référence à l'indice de balayage. [0118] On fait l'hypothèse que les niveaux d'amplitude restent faibles et que l'on cherche principalement à étudier l'influence des micro-impacts sur la stabilité. Alors, on peut supposer que la déformée du mode ainsi que sa pulsation évoluent peu. En d'autres termes, la déformée q,P et la pulsation oi9 n'évoluent pas avec pr . Ainsi, Xr se déduit de : [0119] [Math.27] Xe n lo e : X- 1 [0120] La déformée 13 restant constante, j et, est déduite à partir du calcul des valeurs propres de la matrice A in qui représente la matrice dynamique de l'instabilité au pas d'itération m. [0121] Expression des matrices KI et Cl [0122] Les matrices K t et C t qui représentent respectivement les matrices de raideur et d'amortissement équivalentes sont définies à partir de l'équation Math 14 et peuvent être déterminées analytiquement. En effet, pour chaque élément de contact n, le déplacement relatif normal àn peut s'écrire : [0123] [Math.28] = 6r COSUÙ +6 Sin (w [0124] Par ailleurs, les efforts dans la direction d, qui peut être tangentielle ou normale, entre deux noeuds en contact sont : [0125] [Math.29] [0126] Ainsi, pour chaque direction d, ces efforts doivent vérifier l'équation : [0127] [Math.30] a.rr i t kr' t: ± [0128] La résolution des équations donne : [0129] [Math.31] nal s si= 0 63:r Ôet ' 7.--` 0 ou bs =o sin [0130] Synthèse 101311 En résumé, la démarche générale de la méthode MASA mono-mode peut être syn- thétisée en référence à la figure 2 de la façon suivante. 101321 Etape 200 d' initialisation : on met en oeuvre une détermination de valeurs propres du système de frein par exemple par une méthode CEA classique sur le système initial. Chaque mode instable j est identifié ainsi que sa déformée lie et sa pulsation associée w!) . Cette étape correspond à un facteur d'amplitude modale [0133] Les étapes 210 à 260 sont ensuite mises en oeuvre pour chaque mode instable j identifié. L'indice correspondant à l'itération courante est m. 14 [0134] Etape 210 : on fixe une valeur de l'amplitude modale p71. On effectue un balayage en amplitude modale, c'est-à-dire qu'on considère un vecteur P [P.1 .P , ,p d ,on prend la p-ième valeur du vecteur . Les plages de variations sont à déterminer par l'utilisateur . A partir de l'amplitude modale , représentant la déformée modale non normalisée EP] ;yr est déterminée à partir de l'équation Math 27. 101351 Etape 220: les coefficients de Fourier des déplacements et LI, te 1, ns,m_ sont calculés à partir de l'équation Math 19. aui 101361 Etape 230 : le déplacement tau (t, CULest déterminé en appliquant la transformée de Fourier inverse. 1 [0137] Etape 240 : les efforts non-linéaires fzeit {Au r 0,7 ti),9 )) sont calculés à partir de l'équation Math 7 et exprimés dans l'espace de Fourier. 101381 Etape 250 : les matrices KI et Cr représentant respectivement les matrices de raideur et d'amortissement équivalent de Q sont déterminées à partir de l'équation Math 31. 101391 Etape 260 : la matrice d'état Am pour l'itération m est déterminée et les valeurs propres de cette matrice d'état Am sont déterminées. 101401 Les étapes 220 à 260 sont répétées pour toutes les valeurs d'amplitude modale -.Ir, considérées. [0141] Cas multi-mode (ou multi-fréquence) [0142] Le formalisme de la méthode MASA dans le cas multi-fréquentiel (i.e. multi-mode) va maintenant être présenté. [0143] Dans de nombreux systèmes auto-entretenus, plusieurs instabilités sont identifiées lors d'un calcul par la méthode CEA et les réponses des systèmes dynamiques peuvent être quasi-périodiques. Il est alors nécessaire de développer des méthodes permettant de traiter le cas multi- instable. Cette partie propose une adaptation de la méthode présentée précédemment au cas où Ni instabilités sont identifiées par la méthode CEA. [0144] Écriture des déplacements dans le cas multi-fréquentiel [0145] La solution recherchée ici est supposée quasi-périodique, c'est-à-dire qu'elle est ca- ractérisée par Ni fréquences fi d'excitations incommensurables, c'est-à-dire que les Ni fréquences sont telles que [0146] [Math.32] V i F) [1,Nii2 e N 101471 Dans le cas des systèmes non-linéaires, le contenu fréquentiel des oscillations s'écrit comme une combinaison linéaire de toutes les fréquences fondamentales incommensurables. Ainsi, si Nh harmoniques de la solution sont retenues, alors les différentes fréquences observables sont de la forme : 101481 IMath.331 k1 0}1 (i) [0149] où [0150] [Math.34] c [-Nh, et j e [0151] Ainsi, la réponse du système peut être décomposée comme dans le cas mono- fréquentiel en une série de Fourier généralisée, c'est-à-dire que le déplacement AU peut s'écrire : [0152] [Math.35] [0153] où les coefficients et kl...kNi k k correspondent aux coefficients de Fourier associés à toutes les corn, binaisons linéaires possibles des composantes fréquentielles ml, (02, (Dm Une écriture plus condensée consisterait à introduire la base [0154] [Math.36] [0155] et à poser : [0156] [Math.37] (.0 [0157] Où [0158] [Math.38] tTi T2 44. TNi 101591 est appelée variable d'hyper-temps. Cette variable d'hyper-temps est 2n périodique 16 sur chaque dimension ri. On peut alors définir [0160] [Math.39] [0161] le vecteur du nombre d'harmoniques dans chaque dimension ri. Avec ces notations, l'équation Math 35 peut s'écrire sous la forme : [0162] [Math.40] Ali(t) =U ':..a A kcizN kE. [0163] Il est à noter que cette transformation permet de passer d'un domaine temporel t E 110 un domaine temporel multi-périodique T E RArt où chaque dimension ri est caractérisée par une fréquence fondamentale fi du système dynamique. Ainsi, l'équation Math 40 est 2n périodique dans chaque dimension et AU(r) est une fonction de temps multi-dimensionnelle qui décrit l'évolution temporelle de la fonction multi-fréquentielle AU(t). Cette fonction est souvent appelée fonction tore. Par exemple, dans le cas de la dimension 2, la fonction AU(t) est bi-périodique et évolue dans lle sur un tore formé par les deux paramètres d'hyper-temps (r1, r2). Un exemple est donné en FIG. 3. [0164] L'espace de temps physique est remplacé par un espace multi-dimensionnel dont la paramétrisation suivante permet de trouver r correspondant à un instant t : [0165] [Math.41] t9i TN,, = itSN; [0166] De cc fait, il existe une bijection entre le temps t et les paramètres Dans le cas de la dimension 2, l'illustration du passage du domaine temporel au domaine hyper-temps est représenté en FIG. 4. Ainsi, pour tout instant t, il est possible d'associer un unique point dans l'ensemble [0167] [Math.42] [Ti ATN. [0168] Il n'existe à ce jour aucune méthode d'approximation permettant de faire le passage inverse, c'est-à-dire trouver l'instant t correspondant à un point de l'ensemble défini à l'équation Math 41. La conservation de l'orthogonalité de la série de Fourier permet cependant d'utiliser une transformée de Fourier multi-dimensionnelle dans une 17 procédure AFT. [0169] Ainsi les différentes équations écrites dans le domaine de Fourier mono-fréquenticl restent valables dans le domaine multi-fréquentiel et les déplacements AU et les efforts non- linéaires Q peuvent alors s'écrire sous la forme : [0170] [Math.43] (t) = LZ J0 + c(lcr) [0171] et [0172] IMath.441 flet, [0173] Comme précédemment, seuls les termes de première harmonique sont retenus. C'est-à-dire que : [0174] [Math.45] Ni AU/ = IW kp. [0175] [Math.46] Z-4 iZt nvk C 40 kt) V fiS Sint - t v L [0176] où vk est un vecteur nul de taille N; dontseule la k ième composante est égale 1. [0177] Linéarisation [0178] L'équation du mouvement Math 8 projetée sur les termes de première harmonique s'écrit alors : [0179] IMath.471 [0180] L'orthogonalité de la série de Fourier permet d'écrire les Ni équations suivantes : [0181] [Math.48] VieE it,Nd> +CM:, KnEtt, [0182] et sur chaque harmonique, 5-2,1, est linéarisé de la façon suivante : [0183] [Math.49] il-k. +C AUvk. [0184] Ainsi, on aboutit à Ni sous-systèmes Yk de la forme : [0185] [Math.50] 18 V k N [0186] où [0187] [Math.51] 1: Ativi A tift [0188] et les matrices d'état dynamiques AL sont définies de la façon suivante : [0189] [Math.52] s fui Kyk.) M [0190] Ainsi, à chaque instabilité k est associé un sous-système dynamique sur lequel il est possible de réaliser une analyse de la divergence ou de l'atténuation de la solution dynamique linéarisée autour d'un point d'équilibre pour un jeu d'amplitude modale choisie. L'analyse de la divergence de chaque mode instable k se fait à partir de la partie réelle de la valeur propre associée au mode k dans le sous-système k. En effet, chaque Yk est lié uniquement à (th et il n'y a donc pas de contribution modale des modes j avec j # k dans Yk. De plus, la projection sur la première harmonique ne fait pas apparaitre de terme en Yk. [0191] Dans le cas où un sous-système k ferait émerger une instabilité sur le mode j, pour savoir si la divergence associée au mode j se traduit effectivement par une divergence dans l'équation Math 47 du mouvement, il est nécessaire de considérer le sous-système j qui est le seul à donner la contribution du mode./ dans l'équation du mouvement. Dans le cas où le mode j est identifié stable par la méthode CEA, il est possible de l'ajouter à la base utilisée pour MASA. [0192] Calcul des amplitudes modales [0193] Comme précédemment, on peut écrire AU 1 sous la forme : [0194] [Math.53] [0195] où chaque [0196] [Math.54] 19 [0197] représente une déformée modale non normalisée. [0198] Si les déplacements AU i sont non nuls et connus, alors chaque valeur propre ),..; est déduite de la matrice k associée au sous-système j correspondant de la façon suivante: [0199] [Math.55] IP A [0200] La réponse quasi-périodique d'un tel système est pilotée principalement par les modes instables identifiés linéaiisés autour de la position d'équilibre, c'est-à-dire que l'on peut approcher les déplacements par : [0201] [Math.56] 102021 où les ck sont des coefficients réels. 102031 De la même manière que précédemment, si IN; modes instables de modes propres, correspondant aux déformées : 102041 IMath.571 (Pi 102051 sont identifiés par la méthode CEA avec une valeur propre associée : 102061 IMath.581 [0207] alors on considère les différents vecteurs d'amplitude modale [0208] [Math.59] N [0209] où chaque pel correspond à l'amplitude modale associée au mode k pour le pas de balayage numéro in. On suppose ici aussi que les niveaux d'amplitude restent suffisamment faibles afin que les déformées et les pulsations ny2 n n'évoluent pas avec pro. Ainsi, pour chaque mode k on a : [0210] [Math.60] xen 4 k (fi k>Xne [0211] Puis chaque valeur propre Xi, est obtenue en déterminant les valeurs propres de la matrice Ak. [0212] 1-1 est envisageable dans un premier temps de considérer une grille pour pouvoir effectuer le balayage sur p niais cette approche serait coûteuse numériquement. Pour cette raison, une approche plus restrictive, mais qui propose une vision "physique" du développement des réponses vibratoires est considérée et consiste à revenir dans le temporel afin d'établir un lien entre les différentes amplitudes modales. [0213] Ainsi, si chaque mode est caractérisé par une amplitude modale [0214] [Math.61] [0215] en un instant ln et une valeur propre An alors l'amplitude pk(t',,,) à l'instant +1 tTz t M sera liée à pk(t',1) et à rite de la façon suivante : [0216] [Math.62] fer 1 [0217] Où [0218] [Math.63] [0219] correspond à la fonction partie réelle, cette fonction étant appliquée ici à la valeur propre 102201 Cette équation est équivalente à : 102211 1Math.641 V(i> 102221 Ainsi, les différentes amplitudes modales sont liées les unes aux autres et avec cette approche le balayage ne se fait plus en amplitude modale mais en temps. Le choix du pas de temps est fait de sorte que l'équation Math 62 reste valable, c'est-à-dire que considérer un accroissement constant et égal à : 102231 1Math.651 21 [0224] entre les instants t' et t, est une approximation valide. Il est à noter que chaque valeur propre me- est déterminée à partir de la matrice Ali', associée au sous-système du mode k à l'instant L. [0225] Calcul des efforts non-linéaires Q [0226] Les efforts non-linéaires n sont définis comme précédemment par : [0227] [Math.66] fi .Fni(Us + AU) -- Fni(U,$) - L'AU 102281 L'utilisation d'un algorithme AFT permet de les calculer directement à partir de l'expression de AU(t). Ils sont ensuite projetés sur chaque harmonique et les coefficients de Fourier sont déterminés. 102291 11 est important de remarquer que le couplage entre les différents modes se fait à travers ces efforts non-linéaires. En effet, en considérant une troncature à l'harmonique 1, l'équation Math 66 devient : 102301 IMath.671 1n4 , -Fiti s -.1tzi > Auvk [0231] Et plus précisément, c'est le terme [0232] [Math.68] VIII (US ± k. ..) 102331 qui définit le couplage des différents modes entre eux. 102341 Dans les équations Math 67 et Math 68 F ni représente les efforts non-linéaires de contact et de frottement à l'interface de friction, .1', représente la matrice Jacohienne des efforts non-linéaires autour de la position d'équilibre statique Us représente une position d'équilibre statique, AU'k représente la composante sur la pulsation cuk des déplacements AU. [0235] Synthèse [0236] En résumé, un exemple de réalisation d'un procédé de prédiction d'instabilités de crissement d'un système de frein correspondant à la démarche générale de la méthode MASA multi-mode présentent& en détail ci-dessus peut être synthétisé en référence à la FIG. 5 de la façon suivante. [0237] Bicn que les étapes soient décrites de manière séquentielle, l'homme du métier comprendra que certaines étapes peuvent selon le cas être omises, combinées, exécutées dans un ordre différent et / ou en parallèle. [0238] Etapc 500: on met en oeuvre une méthode d'identification des modes instables au point d'équilibre statique. On utilise par exemple la méthode CEA classique. Les Ni modes instables sont identifiés. La déformée [0239] [Math.69] [0240] et pulsation associée à chaque mode instable est obtenu par cette méthode d'identification de modes instables. [0241] Le vecteur d'amplitude modale défini par Math 59 est initialisé une valeur initiale p° = O. Une valeur de perturbation initiale Ap est choisie. [0242] Un processus itératif est mis en oeuvre. A chaque itération, l'indice m est incrémenté. La valeur de in est fixée à 1, in représentant l'indice du temps ou le numéro de la boucle d'itération dans le processus itératif [0243] Un modèle linéaire du système de frein est généré dans lequel les efforts non- linéaires sont linéarisés vis-à-vis du déplacement et de la vitesse selon des degrés de liberté du système de frein. [0244] En outre, le modèle linéaire du système de frein peut être décomposé en sous- systèmes mécaniques, chaque sous-système correspondant à un mode instable obtenu par la méthode d'identification de modes instables et correspondant à une pulsation ou fréquence. Cela revient, à projeter le modèle linéaire initial sur les fréquences des modes instables identifiés. [0245] Une simulation d'une excitation du système de frein est alors mise en oeuvre en imposant un déplacement particulier selon les degrés de liberté du système de frein. Ce déplacement particulier peut s'exprimer sous forme de combinaison linéaire de N déformées modales, une déformée modale correspondant à un mode instable et étant définie par une déformée modale normalisée, une amplitude modale et une fréquence parmi les N fréquences identifiées par la méthode CEA. [0246] A chaque itération, les niveaux des amplitudes modales varient et sont déterminées en fonction de la variable temps à partir de Math 62. Les déformées permettant l'expression des déplacements imposés ainsi que les fréquences associées qui permettent la mise en oeuvre de la méthode AFT, correspondent aux déformées normalisées et aux fréquences identifiées par la méthode CEA. [0247] Puis on détermine, sur la base du modèle linéaire du système de frein, les variations des valeurs propres des matrices d'état des différents sous-systèmes lorsque les niveaux des amplitudes modales varient. [0248] Une identification d'instabilités est ensuite possible sur la base des valeurs propres complexes ainsi calculées. [0249] Du point de vue mathématique, le processus itératif inclut les étapes 510 à 570. [0250] Etapc 510: l'amplitude modale 23 par est déterminée par récurrence à partir de l'amplitude modale de l'itération précédente pir -1 sur la base de l'équation Math 62 si m>l, ou à partir de la perturbation initiale Ap si tn = 1 [0251] Etape 520 : à partir de chaque amplitude modale Pj , les déformées modales non normalisées X t sont déterminées à partir de l'équation Math 60. [0252] Etape 530 : les déplacements Ativckinet Atjs,k'ni sont calculés à partir des équations Math 53 et Math 54 déclinées pour l'itération m. [0253] Etape 540 : chaque AuFk est déterminé en appliquant la transformée de Fourier inverse aux AUe'm et vk AU 5-m. vk [0254] Etape 550: la somme des efforts non-linéaires le, pour tous les modes k= 1 à Ni est calculée à partir de l'équation Math 67 et exprimés dans l'espace de Fourier. [0255] Etape 560 : les matrices Kos et C.k sont déterminées à partir de l'équation Math 49. [0256] Ce calcul peut s'effectuer de la façon suivante. Pour chaque élément de contact n, la projection sur chaque pulsation wk du déplacement relatif normal s'écrit : [0257] [Math.70] 81: (t) = (5, cos wkt + 8s sin Wkt.. 102581 Par ailleurs, les efforts dans la direction d, qui peut être tangentielle ou normale, entre deux noeuds en contact sont : 102591 IMath.711 nvniid (t) Tc COS (-0k T Ts SiruÉkt 102601 Ainsi pour chaque direction d, ces efforts doivent vérifier l'équation : 102611 IMath.721 fli,r/d(t) = kn"18n(t)+ Cn g et (0 [0262] la résolution des équations donne : [0263] [Math.73] 24 [0264] Etape 570 : les différentes matrices Ank' sont déterminées à partir de l'équation Math 52 et les valeurs propres des différentes matrices .Ait sont déterminées. [0265] Les étapes 510 à 570 sont ensuite répétées pour m= [0266] L'analyse de la divergence de chaque mode instable k se fait à partir de l'analyse de l'évolution en fonction du temps de sa valeur propre 4 dans le sous-système k. Si une divergence du mode .j est observée dans le sous-système k, alors il faut se référer au sous-système j pour voir si cette divergence se traduit effectivement par une divergence dans l'équation du mouvement. Si le mode ,j était identifié comme stable initialement par la méthode CEA, il serait alors nécessaire d'ajouter ce mode à la base MASA afin de pouvoir statuer sur sa stabilité, c'est-à-dire d'ajouter à une itération donnée (m > 1) une fonction de déformation supplémentaire et pulsation associées à ce mode nouvellement identifié comme instable. En alternative, on peut suivre l'évolution temporelle des valeurs propres des N modes instables identifiés au départ par la méthode CEA, puis redémarrer à partir de l'étape 500 une nouvelle simulation avec les N. instabilités supplémentaires qui apparaissent au cours de la simulation afin d'appliquer à nouveau la méthode MASA sur les Ni + Ns modes instables. [0267] Selon le même principe, il serait envisageable de retirer de l'analyse MASA multi- fréquentielle tout mode dont l'amplitude se serait stabilisée et qui donc ne s'avère pas comme étant instable en régime dynamique. On peut alors retirer à une itération donnée (m > 1) la ou les fonctions de déformation et pulsation associées à tout mode considéré comme stable. [0268] La valeur de perturbation Ap en amplitude modale peut être initialisée de différentes manières. Une valeur trop faible implique un plus grand nombre de calculs nécessaires pour arriver à la stabilisation, une valeur trop élevée peut quant à elle conduire à une divergence de la méthode. 102691 Réduction du modèle éléments finis de frein 102701 Afin de pouvoir appliquer la méthode sur un modèle éléments finis de frein in- dustriel, il est possible de réduire la taille du problème. Ainsi, un modèle réduit tel que celui présenté dans le document FAZIO, O., S. NACIVET et J.-J. S1NOU, 2015, « Réduction strategy for a brake svstem with local frictional non-linearities-application for the prédiction of unstable vibration modes », Applied Acoustics, vol. 91, p. 12-24, est utilisé. Dans cette étude, un processus numérique basé sur une réduction modale pour minimiser la taille du modèle est possible grâce à une modélisation non-linéaire spécifique à l'interface de friction. 102711 Comme expliqué dans ce document, la surface de contact plaquette/disque est réduite en un nombre limité de noeuds sur chaque surface de contact. Ces noeuds sont utilisés par la suite en tant que noeud de condensation et en tant que noeud de définition pour des contacts node-to-node. Une fois la surface réduite générée, le modèle réduit est construit en générant un assemblage de Super Elements (SE). Pour une description complète de la génération, de la définition et de l'assemblage des SE, le lecteur pourra se référer aux documents FAZIO, 0.2016, Amélioration de la prédictivité des calculs de crissement de frein, thèse de doctorat, Ecole Centrale de Lyon et FAZIO, 0., S. NACIVET et J.-J. SINOU. 2015, «Reduction strategy for a hrake system with local frictional non-linearities-application for the prediction of unstable vibration modes », Applied Amnisties, vol. 91, p. 12-24 . De plus, une réduction de Craig et Bampton classique est réalisée où les points de réduction de l'interface plaquette/disque sont conservés comme noeuds de réduction. De cette manière, les efforts non-linéaires de contact peuvent être pris en compte dans le présent modèle. [0272] Des exemples de prédiction sont illustrés par référence aux FIGS. 6 à 8 dans des cas dans lesquels un mode instable, non identifiable par la méthode CEA, est détecté grâce à la méthode MASA multi-mode. [0273] Dans l'exemple des FIGS. 6A-6B, la méthode CEA identifie une unique instabilité à 3870 Hz avec une partie réelle de 87. La méthode MASA mono-instable peut être appliquée au départ avec un balayage en amplitude modale. Les évolutions en fonction de l'amplitude modale des parties réelles des valeurs propres et des fréquences obtenues sont représentées respectivement sur la FIG. 6A (parties réelles) et sur la FIG. 6B (fréquences). [0274] On remarque ici que le mode initialement instable (mode MI représenté en trait plein) voit sa partie réelle s'annuler lorsque les amplitudes de vibration augmentent (vers une valeur de 100 d'amplitude modale). Cependant, avant que ce mode MI ne se stabilise, un autre mode M2 instable émerge à 4750 Hz environ (mode représenté en traits pointillés). Ainsi, lorsque l'amplitude modale prend une valeur p = 80, deux modes instables sont présents dans le système et pour pouvoir statuer sur le comportement dynamique global et la contribution des modes instables sur la divergence ou atténuation de la réponse dynamique du système, il est nécessaire d'ajouter ce second mode instable à l'analyse MASA. [0275] Ainsi, une seconde analyse est réalisée avec la méthode MASA multi-mode (avec deux modes dans ce cas précis). La méthode MASA multi-mode est donc initialisée (étape 500) avec le mode Ml instable et le mode M2 qui devient potentiellement instable. Une faible amplitude modale initiale est affectée initialement aux deux modes instables MI et M2. Les évolutions temporelles des amplitudes modales (FIG. 7A), des parties réelles des valeurs propres (FIG. 7B) et des fréquences (FIG. 7C) sont représentées dans les FIGS. 7A-7C. Sur ces figures, les courbes d'évolution des pa- 26 ramètres associés au premier mode instable M1 sont en traits pleins, les paramètres associés au mode instable M2 sont en traits pointillés. [0276] Au départ, seul un mode Ml est instable avec une partie réelle de 87 (FIG. 7B) et une fréquence de 3880 Hz (FIG. 7C). Jusqu'à t = 0.045, sa partie réelle reste sensiblement constante et son amplitude modale augmente. Puis, à partir de t = 0.045 s, le second mode M2 voit la partie réelle dc sa valeur propre augmenter et devenir positive et le premier mode a sa partie réelle qui s'annule rapidement (cf. FIG. 7B). A partir de t = 0.05 s, le premier mode Ml a une partie réelle nulle et s'est stabilisé à une amplitude d'environ 65 (FIG. 7A). Le second mode M2 a sa partie réelle qui reste constante et égale à 21 (FIG .7B) dc t = 0.05 s jusqu'à t = 0.17 s et son amplitude modale augmente progressivement (FIG. 7A). A partir de t = 0.17 s, la partie réelle du second mode M2 augmente à nouveau et croît plus fortement à partir de t = 0.21 s (FIG. 7B). L'amplitude modale associée augmente rapidement (FIG. 7A). Puis à t = 0.22 s, la partie réelle du second mode M2 diminue brusquement et à t= 0.23 s elle s'annule (FIG .7B). L'amplitude modale se stabilise alors à p = 195 (FIG. 7A). La fréquence du second mode M2 est autour de 4730 Hz (FIG. 7C). [0277] Les évolutions des valeurs propres de deux sous-systèmes sont représentées aux FIGS. 8A-8D. Les évolutions des parties réelles des valeurs propres et fréquences correspondantes sont illustrées pour le premier sous-système correspondant au premier mode instable aux FIGS. 8A et 8B et aux FIGS. 8C et 8D pour le second mode instable correspondant au deuxième sous système. [0278] Les évolutions des valeurs propres des deux sous-systèmes de la FIG. 8 montrent qu'aucune autre instabilité potentielle n'émerge dans l'un des sous-systèmes et qu'il n'est donc pas nécessaire d'appliquer la méthode MASA avec un plus grand nombre d'instabilités. [0279] Les FIGS. 8A-8D correspondent aux valeurs propres des matrices AI et A2 des deux sous-systèmes, chaque matrice ayant 4 valeurs propres correspondant à chacun des modes, et on représente l'évolution de ces 4 valeurs propres. La participation d'un mode donné dans la réponse globale est donnée par la valeur propre de ce mode dans le sous-système correspondant. [0280] Pour résumer, sur les FIG. 8A-8D on observe que la méthode MASA indique les contributions suivantes : [0281] t 0.23 : deux modes instables M2 et M4 à 3880 Hz et 4730 Hz à parties réelles nulles qui caractérisent une solution stationnaire quasi-périodique approchée. [0285] Conformément aux résultats présentés sur le modèle académique, il est possible de dire que le système global a une réponse dynamique pilotée par deux modes à 3880 Hz et 4750 Hz avec des amplitudes respectives de 65 et de 195. Ce résultat illustre l'importance de l'utilisation de la méthode MASA multi-modes en complément d'une méthode CEA. En effet, seule une instabilité était détectée avec la méthode CEA alors que deux sont présentes dans la dynamique du système. [0286] Dans certains modes de réalisation, la réponse dynamique stationnaire quasi-pé- riodique du système peut être approchée par : [0287] [Math.75] tiall(t) /hl [0288] avec pl = 65 et p, = 195, = 2.44 104 rads et o = 2.97 104 rad/s et (02 kit et po étant les déformées modales dès deux modes obtenues avec la méthode CEA. 2 [0289] Par ailleurs, ces résultats illustrent aussi l'importance d'utiliser une méthode MASA avec une approche multi-instabilité plutôt que de superposer les résultats de plusieurs approches mono-instables. En effet, dans le cas mono-instable, le premier mode se stabilisait pour une amplitude de 80 alors que l'approche multi-instable prédit une amplitude de 65. [0290] Selon un mode de réalisation, tout ou partie des étapes d'un ou des procédés décrits dans ce document sont mises en oeuvre par un logiciel ou programme d'ordinateur. [0291] Les fonctions et procédés décrits clans ce document peuvent ainsi être mises en oeuvre par logiciel (par exemple, via un logiciel sur un ou plusieurs processeurs, pour exécution sur un ordinateur à usage général (par exemple, via l'exécution par un ou plusieurs processeurs) afin d'implémenter un ordinateur à usage spécifique ou similaire) et/ou puissent être mises en oeuvre dans du matériel (par exemple, en utilisant un ordinateur à usage général, un ou plusieurs circuits intégrés spécifiques (ASIC) et/ou tout autre équivalent matériel). [0292] L'invention concerne ainsi un logiciel ou programme d'ordinateur, susceptible d'être exécuté par un dispositif informatique (par exemple, un ordinateur), au moyen d'un ou plusieurs processeurs de données, ce logiciel / programme comportant des instructions pour cause l'exécution par ce dispositif informatique de tout ou partie des étapes de 28 l'un ou des procédés décrits dans ce document. Ces instructions sont destinées à être stockées dans une mémoire d'un dispositif informatique, chargées puis exécutées par un ou plusieurs processeurs de cc dispositif informatique de sorte à causer l'exécution par ce dispositif informatique du procédé concerné. [0293] Ce logiciel / programme peut être codé au moyen de n'importe quel langage de pro- grammation, et être sous la forme de code source, code objet, ou de code intermédiaire entre code source et code objet, tel que dans une forme partiellement compilée, ou dans n'importe quelle autre forme souhaitable. [0294] Le dispositif informatique peut être mis en oeuvre par une ou plusieurs machines phy- siquement distinctes. Le dispositif informatique peut présenter globalement l'architecture d'un ordinateur, incluant des constituants d'une telle architecture: mémoire(s) de données, processeur(s), bus de communication, interface(s) matérielle(s) pour la connexion de ce dispositif informatique à un réseau ou un autre équipement, interface(s) utilisateur, etc. [0295] Dans un mode de réalisation, tout ou partie des étapes d'un ou plusieurs procédés décrits dans ce document sont mises en oeuvre par un dispositif doté de moyens de mise en oeuvre de ces étapes du ou des procédés concernés. [0296] Ces moyens peuvent comprendre des moyens logiciels (software) (par ex. des ins- tructions d'un ou plusieurs composants d'un programme) et/ou moyens matériels (hardware) (par ex. mémoire(s) de données, processeur(s), bus de communication, interface(s) matérielle(s), etc). [0297] Des moyens mettant en oeuvre une fonction ou un ensemble de fonctions peuvent également correspondre dans ce document à un composant logiciel, à un composant matériel ou bien à un ensemble de composants matériels et/ou logiciels, apte à mettre en oeuvre la fonction ou l'ensemble de fonctions, selon ce qui est décrit ci-dessous pour les moyens concernés. [0298] L'invention concerne aussi un support d'informations lisible par un processeur de données, et comportant des instructions d'un programme tel que mentionné ci-dessus. [0299] Le support d'informations peut être n'importe quel moyen matériel, entité ou dispositif, capable de stocker les instructions d'un programme tel que mentionné ci-dessus. Les supports de stockage de programme utilisables incluent les mémoires ROM ou RAM, les supports de stockage magnétiques tels que des disques magnétiques et des bandes magnétiques, les disques durs ou des supports de stockage de données numériques à lecture optique, etc ou toute combinaison de ces supports. [0300] Dans certains cas, le support de stockage lisible par ordinateur n'est pas transitoire. Dans d'autres cas, le support d'informations peut être un support transitoire (par exemple, une onde porteuse) pour la transmission d'un signal (signal électromagnétique, électrique, radio ou optique) porteur des instructions de programme. Ce 29 signal peut être acheminé via un moyen de transmission approprié, filaire ou non filaire: câble électrique ou optique, liaison radio ou infrarouge, ou par d'autres moyens. Un mode de réalisation concerne également un produit programme d'ordinateur comprenant un support de stockage lisible par ordinateur sur lequel sont stockées des instructions de programme, les instructions de programme étant configurées pour causer la mise en oeuvre par le dispositif informatique de tout ou partie des étapes d'un ou des procédés décrits ici lorsque les instructions de programme sont exécutées par un ou plusieurs processeurs et/ou un ou plusieurs composants matériels programmables. [Revendication I]