Brevet FR2942035A1 — publié le 20100813

PROCEDE DE QUANTIFICATION DE MODES VIBRATOIRES AUTOENTRETENUS. L'invention concerne un procédé pour quantifier des modes vibratoires autoentretenus, notamment de grandeurs physiques d'un système physique non linéaire. Plus particulièrement, l'invention concerne un procédé pour quantifier des modes vibratoires autoentretenus dans un système de frein pour véhicule, notamment pour véhicule automobile. On sait qu'il arrive à des freins de véhicules de crisser à très faible vitesse des roues, par exemple en descendant une rampe d'accès à un parking souterrain. Ce crissement correspond à une fréquence de vibration audible surprenante du fait de l'absence de force d'excitation oscillante. En effet ces vibrations peuvent être observées pour de faibles déplacements angulaires de roue avant même d'avoir atteint une période de révolution. Ce type de vibration ne résulte donc pas d'une excitation périodique. Ce phénomène est attribué à une présence de non linéarités dans le système qui peuvent générer un mode vibratoire autoentretenu tel que le montre par exemple le document de Sinou J-J et al. : � Application of a non linear modal instability approach to brake systems � Journal of Vibration and Acoustics by ASME, janvier 2004, Vol.126, n°1, pages 101 à 107. La méthode d'analyse modale non linéaire complexe qui y est présentée, reste applicable à des systèmes simples avec un nombre de degrés de libertés faible pour lesquels il existe généralement une seule valeur propre à l'origine d'instabilités possibles du système. Cependant, dans un système physique avec de nombreux degrés de libertés, plusieurs modes propres de vibration peuvent coexister avec différentes fréquences incommensurables, c'est-à-dire des fréquences qui ne peuvent pas être considérées comme des harmoniques d'une même fréquence fondamentale. Avec les moyens de traitement de données physiques connus, il devient alors difficile de tenir compte des modes propres, en valeur et en quantité, pour déterminer un mode autoentretenu. L'invention a pour but de quantifier au moins un mode vibratoire autoentretenu de grandeurs physiques d'un système physique non linéaire soumis à un vecteur de force extérieure non oscillante, quelque soit le nombre de modes propres de vibration. Dans ce but, l'invention a pour objet un procédé comprenant : - une étape de réception ordonnancée pour recevoir un mode vibratoire courant comprenant un vecteur de pulsations fondamentales et un vecteur d'amplitudes des grandeurs physiques dans un domaine fréquentiel, pour obtenir un vecteur de déplacement des grandeurs physiques dans un domaine temporel à partir du mode vibratoire courant au moyen d'une transformée inverse de Fourier multifréquence et pour évaluer un vecteur de force non linéaire correspondant audit mode vibratoire courant ; - une étape d'évaluation ordonnancée pour évaluer des valeurs propres correspondant à des instabilités du vecteur de déplacement dans le domaine temporel ; - une étape de production ordonnancée pour produire une fonction du vecteur de pulsations fondamentales, du vecteur d'amplitudes et des valeurs propres évaluées dans le domaine temporel ; et - une étape de test ordonnancée pour calculer une norme de la fonction et pour indiquer que le mode vibratoire courant quantifie le mode vibratoire autoentretenu lorsque ladite norme est inférieure ou égale à une valeur de seuil. Avantageusement, le vecteur de déplacement des grandeurs physiques dans un domaine temporel est obtenu au moyen de la transformée inverse de Fourier multifréquence, sous forme d'un polynôme trigonométrique comprenant des coefficients qui sont les composantes du vecteur d'amplitudes et des monômes trigonométriques dans chacun desquels le temps est multiplié par une combinaison linéaire des pulsations fondamentales. Particulièrement, dans l'étape de production, une matrice élémentaire de raideur dynamique associée à chaque combinaison linéaire des pulsations fondamentales, est remise à jour dans le domaine fréquentiel et les matrices élémentaires de raideur dynamique remises à jour sont regroupées sur une diagonale de matrice globale de façon à obtenir un produit de la matrice globale par le vecteur d'amplitudes dans le domaine fréquentiel. Plus particulièrement, une première partie de la fonction porte sur le produit de la matrice globale par le vecteur d'amplitudes auquel on ajoute le vecteur de force non linéaire et duquel on retranche un vecteur de force extérieure dans le domaine fréquentiel. Plus particulièrement encore, une deuxième partie de la fonction porte sur au moins une partie réelle de valeur propre évaluée dans le domaine temporel. Plus particulièrement aussi, une deuxième partie de la fonction porte sur au moins une différence entre une partie imaginaire de valeur propre évaluée dans le domaine temporel et une pulsation fondamentale. Plus précisément, le procédé comprend une étape d'émission itérative ordonnancée pour réitérer les étapes précédentes en modifiant le mode vibratoire courant de façon à minimiser la fonction tant que ladite norme est supérieure à la valeur de seuil. De préférence, le procédé comprend une étape initiale ordonnancée pour évaluer des valeurs propres correspondant à des instabilités du système physique en équilibre statique, pour définir un mode vibratoire initial comprenant un vecteur initial de pulsations fondamentales en fonction des parties imaginaires des valeurs propres et un vecteur initial d'amplitudes des grandeurs physiques en fonction de vecteurs propres associés aux valeurs propres. Particulièrement, le système physique est un système de frein pour véhicule. L'invention a aussi pour objet un produit programme d'ordinateur chargeable directement dans la mémoire interne d'un calculateur numérique, comprenant des portions de code logiciel pour l'exécution des étapes du procédé selon l'invention lorsque le programme est exécuté sur un calculateur. La description explicative qui va suivre fait référence aux dessins schématiques annexés donnés uniquement à titre d'exemple et illustrant plusieurs modes de réalisation de l'invention et dans lesquels : la figure 1 est un schéma de maillage pour décrire formellement par éléments finis des composants de frein ; - la figure 2 est un organigramme comprenant des étapes d'un mode de mise en oeuvre du procédé selon l'invention pour quantifier des amplitudes et des pulsations de vibrations autoentretenues ; - la figure 3 est un organigramme comprenant des étapes d'un autre mode de mise en oeuvre du procédé selon l'invention pour quantifier des amplitudes et des pulsations de vibrations autoentretenues. La figure 1 montre un exemple de modèle de système mécanique pour véhicule, notamment pour véhicule automobile, auquel l'invention s'applique. Le système mécanique illustré comprend un disque de frein 1 et une plaquette de frein 2. La plaquette 2 est destinée de manière connue à être appliquée contre une face du disque avec une pression suffisante pour générer un frottement entre les faces en contact du disque 1 et de la plaquette 2. Ce frottement assure alors un freinage par transformation d'énergie cinétique en travail des forces de frottement et dissipation de ce travail sous forme de chaleur. Le système mécanique illustré est simplifié et a uniquement pour but d'expliquer comment l'invention est mise en oeuvre. Une mise en oeuvre réelle est plus complexe mais développe les mêmes moyens que ceux ici expliqués. Ainsi, on sait qu'un système de frein comprend généralement une deuxième plaquette, non représentée et destinée à être appliquée symétriquement par rapport à la première plaquette 2 contre la face opposée du disque 1. Le système de frein comprend encore d'autres accessoires non représentés tels qu'un piston hydraulique combiné avec un étrier pour appuyer sur les plaquettes en phase de freinage. Le modèle simplifié représenté sur la figure 1 est du type de modèles plus complets généralement produit par des équipements de CAO (conception assistée par ordinateur) qui sont utilisés pour fabriquer les pièces du système mécanique considéré. Pour étudier le comportement mécanique du système mécanique, il est habituel de recourir à un procédé de simulation par éléments finis. Chaque pièce mécanique est virtuellement découpée en une multitude finie de volumes élémentaires qui constituent chacun un élément fini 10. Chaque élément fini 10 est par exemple un tétraèdre dont chacun des quatre sommets 11, 12, 13, 14 constitue un noeud de maillage partagé avec un ou plusieurs autres éléments finis. Le comportement de chaque pièce résulte d'une déformation de ses éléments finis sous l'effet de contraintes internes provoquées par les éléments finis voisins et de contraintes externes provoquées par des forces extérieures entre les pièces éléments finis se noeuds 11-17,..., tridimensionnel à 3N degrés de associés à un vecteur U (t) composantes, associés chacun distinct. au système et des forces de contact du système. Les déformations des traduisent par des déplacements des N 21-27,... du maillage qui, dans l'espace qui nous est familier, donnent naissance libertés. Les 3N degrés de liberté sont comprenant N triplets de à un déplacement de noeud Une mise en contact de la face inférieure de la plaquette 2 avec la face supérieure du disque 1, provoque une superposition de noeuds 25, 26, 27, de la plaquette 2 avec des noeuds 15, 16, 17, du disque 1. Lorsque la plaquette 2 est appuyée contre le disque 1 par le piston (non représenté), les noeuds 15-17 et 26-27 subissent des forces de contact Fcontact sensiblement perpendiculaires aux faces en contact et des forces de frottement Ffriction sensiblement tangentes à la surface et orientées dans le sens inverse au sens de rotation du disque 1. Pour chaque couple de noeuds 15, 25 du disque et de la plaquette en contact, les forces de contact et de frottement sont exprimées par les relations : F ontact =0 sauf z ù u15 z > 0) = Fcontact, z = a(u25, ù u15, z)+ (3(u Ffriction,x = Fcontact signe(i125,x ùu15,x) Ffricton,y = Fcontact signe(u25,y ùu15,y) où a et R sont respectivement un coefficient de dureté linéaire et un coefficient de dureté non linéaire, où p est un coefficient de frottement, où u25,z et u15,Z sont respectivement une composante du déplacement du noeud 25 et une composante du déplacement du noeud 15 perpendiculairement aux faces en contact, et où Û25,X , Û25,y, Û15,x et ü15,y sont respectivement les composantes de la vitesse de déplacement du noeud 25 et les composantes de la vitesse de déplacement du noeud 15 dans le plan des faces en contact. On observe que les noeuds des faces en contact sont soumis à des forces non linéaires qui sont regroupées dans un vecteur Fn1(U) à 3N composantes dont seules les composantes associées à des composantes de noeuds appartenant à des faces en contact du disque 1 et de la plaquette 2, sont non nulles. D'autres noeuds du maillage sont soumis à des forces 35 extérieures de nature linéaire. Ce sont par exemple les l3 z ùU15,z) noeuds 21, 22 situés sur la face de la plaquette 2 qui est opposée à la face en contact avec le disque 1. Ce sont par exemple encore les noeuds 18, 19 du disque qui subissent le couple d'une roue (non représentée) sur laquelle le disque est fixé. Les forces extérieures sont regroupées dans un vecteur Fext à 3N composantes parmi lesquelles les composantes non nulles sont celles associées à des noeuds qui subissent les forces extérieures. Le comportement temporel du vecteur U(t) à 3N composantes de déplacement, est régi par les lois de la dynamique qui sont exprimées par l'équation suivante : MU(t) + CU(t) + KU(t) + Fä (U) = Fext où M est une matrice de masse à 3N lignes et 3N colonnes, 0(t) est une dérivée seconde du vecteur U(t) par rapport au temps qui comprend les accélérations des noeuds, C est une matrice d'amortissement à 3N lignes et 3N colonnes, U(t) est une dérivée première du vecteur U(t) par rapport au temps qui comprend les vitesses des noeuds et K est une matrice de raideur à 3N lignes et 3N colonnes. Dans la formule ci-dessus énoncée, le vecteur FnI(U) induit des solutions vibratoires auto entretenues pour le vecteur U(t). En d'autres termes, on peut observer des déplacements oscillants de noeuds sans qu'il existe de forces oscillantes d'excitation. Les vibrations engendrées par les forces non linéaires peuvent atteindre des fréquences et des amplitudes audibles, souvent de nature désagréable. Le but du procédé selon l'invention est de prédire, dès la conception, des valeurs de fréquences et d'amplitudes qui pourront se produire de façon à pouvoir juger des modifications à apporter au système avant de le lancer en fabrication. L'invention concerne un procédé pour produire, à partir d'une description formelle d'un système physique, des modes de vibration dynamiquement stables, en d'autres termes autoentretenus, sur lesquels le système physique risque d'accrocher. Une description formelle est une structure matérielle de données qui décrit des propriétés physiques d'un système physique de façon à pouvoir être traitée par une machine, c'est-à-dire de façon automatique sans recourir à un exercice d'activité intellectuelle pour effectuer les traitements qui, de par leur volume, dépassent les capacités cérébrales humaines. L'intelligence humaine intervient par contre sur la conception du procédé et de la machine. Pour des traitements de données volumineux, la machine est souvent réalisée au moyen d'un ordinateur et de programmes informatiques exécutables par l'ordinateur. La figure 1 nous a permis d'expliquer ce qu'est une description formelle de système physique au moyen d'un maillage par éléments finis. Les éléments finis constituent un exemple de description formelle particulièrement adaptée aux systèmes mécaniques tels qu'en l'occurrence un système de frein. L'invention est capable de traiter d'autres descriptions formelles. Les systèmes électriques, par exemple, peuvent comprendre des composants non linéaires tels que des diodes ou des portes logiques à effet de seuil. On connaît des descriptions formelles nommées listes de connectique (netlist en anglais) utilisées pour fabriquer des circuits électriques et auxquelles l'invention est applicable. L'équation de la dynamique ci-dessus mentionnée est alors remplacée par une équation du type de celle dans laquelle la grandeur physique qui est le vecteur déplacement est remplacé par la grandeur physique qui est un vecteur de charges électriques dont les dérivées premières sont des courants et les dérivées secondes des dérivées de courant par rapport au temps, la matrice de masse étant alors remplacée par une matrice d'inductance, la matrice d'amortissement par une matrice de résistance et la matrice de raideur par une matrice de capacitance, les vecteurs de forces par des vecteurs de tensions électriques. La figure 2 est un ordinogramme pour expliquer un premier mode de mise en oeuvre de l'invention. Une étape d'initialisation 100 est ordonnancée pour déterminer des valeurs de départ d'amplitudes et de pulsations de modes vibratoires à partir desquelles le procédé va chercher à converger vers les amplitudes et les pulsations de modes vibratoires autoentretenus qui sont solution de l'équation exprimant les lois de la dynamique. Les valeurs de départ d'amplitudes et de pulsations de modes vibratoires peuvent être choisies de différentes manières. De préférence, elles sont basées sur les valeurs propres représentatives d'un équilibre statique instable du système, en tenant compte des forces non linéaires. Cette heuristique repose sur une hypothèse physique des inventeurs selon laquelle les modes vibratoires en équilibre dynamique sont conditionnés par les équilibres statiques instables. On rappelle que si un équilibre statique correspond à un vecteur déplacement de valeur constante, un équilibre dynamique correspond à un vecteur déplacement oscillant autour d'une valeur constante dans une plage d'amplitudes bornée. Un déplacement Uéqu correspond à une position d'équilibre s'il vérifie, sous forme d'équation d'équilibre, l'équation exprimant les lois de la dynamique avec des dérivées première et seconde par rapport au temps nulles : KUéqä(t)+F1(Uéqä) =Fext L'équilibre est instable lorsqu'une perturbation infinitésimale a tendance à éloigner définitivement le système de son équilibre statique. Il convient donc de mesurer le comportement de Uéqu+dU(t) dans l'équation comportementale: MUéqu + dU(t) + CUéqu + dU(t) + K[Uéqu + dUkt) + F,,, (Uéqu + dU, d U) = Fext En retranchant l'équation d'équilibre de l'équation comportementale, on obtient l'équation des perturbations: MdU(t) + CdU(t) + KdU(t) + Fn1(Uéqu + dU, d U) ù Fn1(Uéqu) = 0 L'équation des perturbations peut encore s'écrire sous forme d'équation d'état : I dU 0 + - m°-'K - M-1C dU_ F l (Uéqu + dU, du) ù Fi (Uéqu ) Où I est la matrice identité. Pour des perturbations {dU,dÛ}'' de valeur faible : 1:n1 (Uéqu +dU,dU dU (Uéqu +Jnl(Uéqu) Fi ) d u' Où JnI (Uéqu) est la matrice jacobienne des forces non linéaires, considérée sur l'état d'équilibre statique. Considérant les vecteurs d'état et de perturbation d'état . U dU X= . et dX= . U dU Et définissant une matrice de gain A telle que : 0 I A= -M-1K -M-1C L'équation des perturbations a la forme suivante au premier ordre . dX = [A+JnI(Uéqu) = J0 (Uéqu )dX dU dU Les valeurs propres A(Jo) de la matrice Jo sont des nombres complexes dont la représentative d'une fréquence est représentative d'un partie imaginaire est et dont la partie réelle équilibre statique asymptotiquement stable lorsqu'elle est négative et d'un équilibre statique instable lorsqu'elle est positive. Le procédé comprend dans l'étape 100, les opérations qui consistent à générer la matrice A à partir de données fournies par un équipement de CAO, à déterminer la matrice jacobienne Jn1 (Uéqu) des forces non linéaires pour les valeurs d'équilibre du système examiné et à calculer les m valeurs propres À(Jo) de la matrice Jo égale à la somme des matrices A et Jni (Uéqu) Les opérations qui suivent dans l'étape 100, sont celles qui consistent à retenir parmi les m valeurs propres À(Jo), les n valeurs propres Ai(Jo) dont la partie réelle est positive ou nulle de façon à élaborer une liste des pulsations Ooi qui constituent chacune la partie imaginaire distincte d'une valeur propre A~(Jo) dont la partie réelle est positive ou nulle. Dans leur approche heuristique pour faire traiter par une machine, les phénomènes vibratoires autoentretenus qui sont provoqués par des non linéarités de système physique, les inventeurs ont élaboré un moyen de traitement qui consiste à exprimer le vecteur déplacement U(t) sous forme d'une série limitée de fonctions oscillantes comprenant un nombre n de pulsations de base w avec leurs premières harmoniques de pulsations kiwi jusqu'à un ordre Nh. La valeur de Nh est choisie dans une plage de valeurs de limite inférieure suffisamment haute pour reproduire correctement le phénomène tel qu'il est observé dans le domaine technique considéré et de limite supérieure suffisamment basse pour ne pas faire exploser les temps de calcul en fonction de l'ordinateur choisi. Par exemple à titre illustratif et non limitatif . U(t)= L Uk cos(klw )t + Uksin(klw )t kE Z =1 Où désigne le produit scalaire d'un vecteur k de n composantes ki par un vecteur i de n composantes coi. On définit alors le vecteur harmonique Û qui comprend les 2N,+1 composantes Uk et Uk pour les N, combinaisons possibles de nombres entiers relatifs ki dans k telles que dans cet exemple: 1 Nh D'autres combinaisons de nombres entiers relatifs ki sont possibles telles que par exemple celles qui consiste à combiner l'ensemble des valeurs entières comprises entre - Nh et + Nh. Nous voyons que dans cet autre exemple, le nombre de termes N, = 2Nh+l! (factorielle de 2Nh+l) est rapidement nettement supérieur à la valeur du premier exemple. En référençant chaque combinaison de ki par un indice q variant de 0 à N,, associé au vecteur k correspondant de sorte que l'indice q=0 correspond au vecteur k dont toutes les composantes ki sont nulles, la série limitée peut encore s'écrire sous forme d'un polynôme trigonométrique: N, U(t)=Uô +LUgcos(k;w )t+Ugsin(kiw )t q=1 Les coefficients de ce polynôme trigonométrique sont les composantes du vecteur harmonique ou d'amplitudes Û. Dans chacun des monômes trigonométriques, le temps est multiplié par une combinaison linéaire des pulsations fondamentales co à coefficients entiers relatifs (integer en anglais). On définit ainsi une transformée de Fourier multiple ou multifréquence FFTM telle que Û = FFTM,, (U (t)) . A la différence de la transformée de Fourier rapide FFT classique dans laquelle chaque pulsation harmonique est un multiple entier d'une unique pulsation fondamentale, les pseudo pulsations harmoniques de la FFTM sont incommensurables au sens où le rapport entre deux pulsations pseudo harmoniques n'est pas nécessairement ni un nombre entier, ni même un nombre rationnel. A la différence de la transformée de Fourier rapide classique qui donne une fonction périodique, la transformée de Fourier multifréquence définie ci-dessus peut donner une fonction pseudopériodique. Dans l'étape 100, le vecteur i est initialisé au vecteur coo de composantes cOoj ci-dessus déterminées. Le vecteur harmonique Û est initialisé au vecteur harmonique Ûo dont les composantes U, et Uq sont nulles sauf celles dont l'indice q est associé à un vecteur k unitaire dont une composante ki et une seule, est égale à 1. Chaque couple de composantes {u ,US}T du vecteur harmonique Ûo est généré au moyen de la formule suivante : US :=pel'i+`P.) avec je [1,n] Où p est une constante quelconque déterminée a priori. Les vecteurs `Ifs et `Ifs sont le vecteur propre et le vecteur propre conjugué de la matrice Jo correspondant à la valeur propre Ai = ri + ici . Après initialisation du vecteur multi fréquences multi harmoniques Û au vecteur Ûo et des pulsations de base co aux pulsations (Ooj dans l'étape 100, une séquence d'étapes 101 à 106 est réitérée jusqu'à obtenir un vecteur multi fréquences multi harmoniques solution Ueolution et des pulsations de base cosolution,j qui donnent les amplitudes et les fréquences d'oscillations autoentretenues potentielles qui sont alors éditées dans une étape 107. L'étape 101 est ordonnancée pour recevoir un mode vibratoire courant qui comprend le vecteur de pulsations fondamentales et le vecteur d'amplitudes des grandeurs physiques dans le domaine fréquentiel tels qu'ils sont émis à partir de l'étape 100 lors d'une première exécution de l'étape 101 et tels qu'ils sont émis à partir d'une étape 106 décrite ci-après lors d'exécutions suivantes de l'étape 101. Du code logiciel permet de calculer le vecteur FnI(U) des forces non linéaires dans le domaine temporel et le vecteur Fnl(Û) des forces non linéaires dans le domaine fréquentiel. Le vecteur U(t) dans le domaine temporel est obtenu par transformée de Fourier multiple inverse appliquée au vecteur Û du domaine fréquentiel. Le vecteur FnI(U) des forces non linéaires est calculé dans le domaine temporel. Dans le cas où le système physique est un système mécanique de freinage, les forces non linéaires comprennent les forces de contact et de frottement exprimées pour chaque couple de noeuds 15, 25 du disque et de la plaquette en contact par les relations : F ntact =0 sauf ù z u15,z > 0) Fcontact,z = a(u z ù u15,z) + z (3(u ù u15,z )3 25, 25, Ffriction,x = Fcontact signe(ù25,x -u15,x) Ffriction,y = Fcontact signe(U25,y -u15,y) Où u15 et u25 sont des triplets de composantes du vecteur U(t) obtenu par transformée de Fourrier multiple inverse dans l'itération courante. Le vecteur des forces non linéaires Fnl(Û) dans le domaine fréquentiel est ensuite obtenu par transformée de Fourrier multiple du vecteur Fn1 (U) . L'étape 102 est ordonnancée pour calculer les valeurs propres À(J) qui résultent des forces non linéaires Fnl (U) qui sont générées par le vecteur déplacement U(t) considéré dans l'étape 101. La matrice Jn1 est régénérée à chaque réitération de l'étape 102 en recalculant les dérivées partielles des forces non linéaires par rapport aux déplacements, prises aux valeurs de déplacement qui sont données par le vecteur U(t) obtenu dans l'étape 101. La matrice A qui, ne dépendant que des masses, des raideurs et des coefficients d'amortissement, ne dépend pas des déplacements, reste de valeur constante à chaque réitération de l'étape 102. La matrice A est donc celle qui est déterminée dans l'étape 100 d'initialisation. La matrice Jn1 et additionnée à la matrice A pour obtenir la matrice J. Les valeurs propres À de la matrice J sont alors calculées au moyen d'une méthode classique de calcul de valeurs propres. De façon à évaluer les valeurs propres qui correspondent à des instabilités du vecteur de déplacement tel qu'il est considéré dans le domaine temporel lors de l'exécution courante de l'étape 101, parmi les valeurs propres courantes qui viennent d'être calculées dans l'étape 102, seules sont retenues celles qui sont les plus proches des valeurs propres retenues précédemment. Ainsi, lors d'une première exécution de l'étape 102, les valeurs propres retenues sont celles qui sont les plus proches des valeurs propres retenues dans l'étape d'initialisation 100, c'est à dire aux valeurs propres à partie réelle positive ou nulle. Lors des exécutions successives de l'étape 102, les valeurs propres retenues calculées, sont celles valeurs propres retenues l'étape 102. L'étape 103 est fonction objective f parmi celles qui viennent d'être qui sont les plus proches des dans une exécution précédente de ordonnancée pour produire une transmise a un outil 108 d'optimisation utilisé pour fournir de nouvelles valeurs de vecteur de déplacement Û et de pulsations cotant que la solution n'est pas trouvée. La fonction objective f comprend une concaténation de deux sous fonctions objectives fi et f2. Dans le mode de mise en oeuvre illustré par la figure 1, la sous fonction objective f2 comprend les parties réelles des valeurs propres retenues dans l'étape 102. La sous fonction objective fi quant à elle, est indépendante de l'étape 102. Alors que l'étape 102 est orientée vers une appréciation du degré d'instabilité d'un équilibre représenté par le vecteur Û qui est considéré en étape 101, la sous fonction objective fi est orientée vers une appréciation de la réponse en fréquence du système physique considéré. Reprenant l'équation qui exprime les lois de la dynamique dans le domaine temporel : MU(t) + CÙ(t) + KU(t) + Fä1 (U) = Fext Reprenant l'expression du vecteur déplacement sous 15 forme d'un polynôme trigonométrique : N, U(t) = Uô + L Ugcos(kiw )t + Ugsin kiwi )t q=1 Les lois de la dynamique s'expriment encore de la façon suivante : KUô + E (- kiwi ) M + K)Uq + kiwi )CUg 4s(kiwi)t q=1 +L(( kiwi)2M+K)Uqù (kiwi )CUqin(kiw )t+F1(t,U,U)ùFext(t)=0 q=1 20 Pour chacune des Ne combinaisons de composantes ki du vecteur k considérée, ou en d'autres termes pour chacune des combinaisons d'harmoniques loci, la machine conforme à l'invention remet à jour une matrice de raideur dynamique Aq associée à la gième combinaison d'harmoniques, en 25 fonction des pulsations co considérées dans l'étape 101 : 2 -(kiwi) M + K (k~wi)C 2 -(kiwi)C M + K Les matrices Aq pour q variant de 1 à Ne, sont regroupées dans une matrice A(co) telle que : Aq :_ K 0 0 0 0 Al 0 0 0 0 0 0 0 0 A A(w;)n 1 := Nc La machine calcule la transformée de Fourier multifréquence Fxt des forces extérieures. Généralement, les forces extérieures étant constantes, la transformée de Fourier multifréquence Fxt des forces extérieures peut être calculée une fois pour toutes dans l'étape 100. La machine utilise ensuite le vecteur des forces non linéaires FI(U) déterminé dans l'étape 101 pour le domaine fréquentiel, de façon à instancier la fonction f1: f (U, ci , J E [1, n ]) := A(cw1)n U+ (U) ù F t On remarque que le vecteur Û de déplacement est 15 conforme aux lois de la dynamique lorsque : f (U,wj,J E [1,n]):= 0 Chercher une solution Û en calculant la fonction fi à chaque itération de l'étape 103 jusqu'à annuler la fonction f1, ressemble à la méthode de la balance 20 harmonique appliquée à un système excité par une oscillation extérieure de pulsation G)ex. On constate cependant plusieurs différences d'importance fondamentale. La méthode de la balance harmonique telle qu'elle est connue est applicable avec une pulsation 25 unique 0ex alors que plusieurs pulsations co sont traitables dans le cadre de l'invention. Dans la méthode de la balance harmonique connue, la valeur de la pulsation wex est identique à chaque itération alors que dans le cadre de l'invention, la valeur des pulsations co varie à chaque itération en fonction des valeurs considérées dans l'étape 101. La machine transmet la fonction f à un outil d'optimisation 108 en charge de trouver un couple de valeurs . qui rapproche la fonction f de zéro. Le couple solution est le couple qui, mettant la fonction fi à zéro, satisfait les lois de la dynamique et qui, mettant la fonction f2 à zéro, est représentatif d'un équilibre dynamique. Ce problème est un cas particulier d'un problème mathématique plus vaste qui est l'optimisation de fonctions objectives. Très peu de systèmes d'équations non linéaires peuvent être résolus analytiquement. Il faut donc obligatoirement recourir à des méthodes numériques itératives. Différentes méthodes itératives sont connues de l'homme du métier. Pour faciliter la compréhension de l'exposé, éléments suivants. De nombreuses méthodes dérivent de la méthode de d'un point initial sur la équations on se déplace de minimum local de la fonction, rappelle simplement les d'optimisation locale à partir surface définie par les manière logique vers un c'est-à-dire en suivant la on Newton dans laquelle, pente de la fonction, comme un skieur qui descendrait la piste vers le télésiège. Tout les pistes en ligne droite comme le skieur qui descend ou avec des virages, les algorithmes d'optimisation ont les mêmes possibilités. Les algorithmes Hybrid et Levenberg-Marquardt semblent être les plus robustes et les plus utilisés dans les logiciels disponibles dans le commerce (MatlabTM, ScilabTM) ainsi que dans les librairies mathématiques gratuites sur Internet (GNU Scientific Library, Minpack, Levmar). L'outil d'optimisation 108 étant connu par ailleurs, il n'est pas nécessaire d'en décrire le fonctionnement en détails. Les explications données à propos de l'outil 108 servent simplement à comprendre quelles en sont les fonctionnalités utilisées, pour faciliter la compréhension de l'invention. L'étape 104 est ordonnancée dans l'outil d'optimisation 108 pour comparer une norme de la fonction f à une valeur de seuil faible c. Différents types de normes peuvent être utilisées telles que la plus grande valeur absolue des composantes de la fonction f, la racine carrée de la somme des carrés des composantes ou encore bien d'autres normes connues. La valeur faible E est prédéterminée de manière connue dans le domaine en tenant compte de la précision du résultat souhaitée, de la précision disponible des calculs et du temps de calcul accepté pour constituer une tolérance sur l'atteinte du nombre zéro. Si la norme de la fonction f est inférieure à E, une étape finale 107 est activée, sinon une étape 105 est activée. Dans l'étape 107, la solution est annoncée comme étant celle qui comprend le vecteur Û et l'ensemble des pulsations co qui ont été utilisés pour la dernière exécution de l'étape 103. Dans l'étape 105, un gradient de descente s(Û,wi) est calculé à partir de dérivées de la fonction f par rapport à ses variables. Dans l'étape 106, le gradient de descente s (Û, wi) est multiplié par un pas de convergence y puis additionné à la valeur précédente du couple {0,wil de façon à obtenir une nouvelle valeur du couple {0,w il avec laquelle les étapes 101 à 103 sont réitérées. Le pas de convergence y est généralement basé sur des considérations balistiques. De manière connue, le pas de convergence y est élevé en début de processus de façon à se rapprocher rapidement de la solution puis diminue progressivement au fur et à mesure que l'on se rapproche de la solution de façon à ne pas passer au-delà. La figure 3 permet d'expliquer un autre mode de mise 5 en oeuvre de l'invention. L'étape d'initialisation 200 est semblable à l'étape 100 mais différentes stratégies d'initialisation sont disponible de façon à pouvoir réitérer l'étape d'initialisation avec d'autres valeurs si nécessaire. Par 10 exemple plusieurs valeur p de multiplication des vecteurs propres sont disponibles ou la somme des vecteurs propres conjugués est remplaçable par leur différence. Les étapes 201 et 202 sont identiques aux étapes 101 et 102. 15 L'étape 203 est comparable à l'étape 103 mais elle en diffère en ce que la fonction f2 ne donne plus les partie réelles des valeurs propres retenues mais les différences entre les parties imaginaires des valeurs propres À retenues et les pulsations code même indice j. 20 Cette différence n'est pas nulle a priori car la partie imaginaire d'une valeur propre À n'est pas nécessairement égale à la pulsation code même indice j à chaque itération de l'étape 203. Il est vrai que lorsque la stratégie d'initialisation expliquée en référence à la 25 figure 2 est utilisée dans l'étape 200, la partie imaginaire de chacune des valeurs propre À est égale à la pulsation coi dans l'étape 200 par construction. Cependant dans l'étape 203, les valeurs propres sont celles de la matrice J déterminée dans l'étape 202 alors que les 30 pulsations coi sont celles considérées dans l'étape 201. Bien que les théories connues des cycles limites et des bifurcations de Hopf à propos des systèmes auto excités ou autonomes enseignent que c'est la présence d'au moins une partie réelle nulle de valeur propre qui 35 est nécessaire et suffisante pour observer des oscillations autoentretenues, les inventeurs ont eu l'intuition qu'une égalité entre partie imaginaire de valeur propre de la matrice J et pulsation d'oscillation du système pouvait être aussi une condition nécessaire et suffisante pour observer des oscillations autoentretenues. L'outil 209 d'optimisation est semblable à l'outil 108. Les étapes 204 à 206 sont semblables voire identiques aux étapes 104 à 106. En effet ce sur quoi porte la fonction f est totalement transparent pour l'outil 108 ou 209. Cependant, en absence de théorie scientifique connue pour corroborer une telle démarche, les inventeurs ont préféré ajouter une étape 207 par sécurité, déclenchée simplement en fin de processus lorsque le test sur la valeur infinitésimale de la norme de la fonction f(Û,wi) est positif dans l'étape 204. Considérant un vecteur constitué des parties réelles des valeurs propres déterminées lors de la dernière exécution de l'étape 202, l'étape 207 teste si la norme de ce vecteur tend vers zéro de façon à vérifier que la solution trouvée est conforme au critère scientifique d'équilibre dynamique oscillant. Une étape 208 identique à l'étape 107 est activée pour éditer la solution fournie en sortie de l'outil 209 en cas de réponse positive au test de l'étape 207. Un rebouclage sur l'étape 200 est prévu pour réitérer l'ensemble du processus avec d'autres conditions initiales en cas de réponse négative au test de l'étape 207. Sur les nombreux essais pratiqués, les inventeurs ont toujours constaté une réponse positive dès le premier test de l'étape 207. Les valeurs d'amplitudes et de pulsations des oscillations autoentretenues qui sont obtenues au moyen du procédé précédemment décrit, sont alors comparées à des seuils de tolérance admissibles, notamment pour l'oreille humaine dans le cas du crissement des freins. En cas de dépassement des seuils, on modifie certains paramètres du modèle auquel on applique alors à nouveau le procédé de l'invention jusqu'à vérifier que les valeurs quantifiées sont ramenées dans des limites acceptables. Dans le cas des freins, par exemple, on modifie le coefficient de frottement p par traitement de la surface du disque ou des plaquettes, on modifie la géométrie ou on ajoute des éléments d'amortissement. Les modifications sont faites au niveau de la modélisation du système physique sans avoir à fabriquer de prototypes pour effectuer des essais, ce dont il résulte des économies de temps, d'énergie et de matière. Lorsque les modifications aboutissent à des valeurs concluantes telles qu'elles sont quantifiées par le procédé, on peut alors lancer le système physique en fabrication de sorte que le système physique produit ne vibrera que dans des tolérances admissibles. Les explications qui viennent d'être données font souvent référence à la pulsation d'une oscillation. La pulsation est une grandeur équivalente à la fréquence au coefficient 2n près. Les explications restent donc valables lorsque l'on remplace le terme pulsation par le terme fréquence implicitement multiplié par 2n. REVENDICATIONS1. Procédé pour quantifier au moins un mode vibratoire autoentretenu de grandeurs physiques d'un système physique non linéaire soumis à un vecteur de force extérieure non oscillante, comprenant : - une étape (101) de réception ordonnancée pour recevoir un mode vibratoire courant comprenant un vecteur de pulsations fondamentales (c~~) et un vecteur d'amplitudes (Û) des grandeurs physiques dans un domaine fréquentiel, pour obtenir un vecteur de déplacement (U(t)) des grandeurs physiques dans un domaine temporel à partir du mode vibratoire courant au moyen d'une transformée inverse de Fourier multifréquence et pour évaluer un vecteur de force non linéaire (Fnl(U), Fnl(Û)) correspondant audit mode vibratoire courant ; - une étape (102) d'évaluation ordonnancée pour évaluer des valeurs propres (À) correspondant à des instabilités du vecteur de déplacement (U(t)) dans le domaine temporel ; - une étape (103) de production ordonnancée pour produire une fonction (f) du vecteur de pulsations fondamentales (co), du vecteur d'amplitudes (Û) et des valeurs propres (À) évaluées dans le domaine temporel ; et - une étape (104) de test ordonnancée pour calculer une norme de la fonction (f) et pour indiquer que le mode vibratoire courant quantifie le mode vibratoire autoentretenu lorsque ladite norme est inférieure ou égale à une valeur de seuil (s). 2. Procédé selon la revendication 1, caractérisé en ce que le vecteur de déplacement (U(t)) des grandeurs physiques dans un domaine temporel est obtenu au moyen de la transformée inverse de Fourier multifréquence, sous forme d'un polynôme trigonométrique comprenant des coefficients qui sont les composantes du vecteurd'amplitudes (Û) et des monômes trigonométriques dans chacun desquels le temps est multiplié par une combinaison linéaire des pulsations fondamentales ((oj). 3. Procédé selon la revendication 2, caractérisé en ce que dans l'étape (103) de production, une matrice élémentaire de raideur dynamique associée à chaque combinaison linéaire des pulsations fondamentales (GO, est remise à jour dans le domaine fréquentiel et les matrices élémentaires de raideur dynamique remises à jour sont regroupées sur une diagonale de matrice globale (A (coi)) de façon à obtenir un produit de la matrice globale (A (coi)) par le vecteur d'amplitudes (Û) dans le domaine fréquentiel. 4. Procédé selon la revendication 3, caractérisé en ce qu'une première partie (f1) de la fonction (f) porte sur le produit de la matrice globale (A (coi)) par le vecteur d'amplitudes (Û) auquel on ajoute le vecteur de force non linéaire Fnl(Û) et duquel on retranche un vecteur de force extérieure Fext dans le domaine fréquentiel. 5. Procédé selon la revendication 4, caractérisé en ce qu'une deuxième partie (f2) de la fonction (f) porte sur au moins une partie réelle de valeur propre (À) évaluée dans le domaine temporel. 6. Procédé selon la revendication 4, caractérisé en ce qu'une deuxième partie (f2) de la fonction (f) porte sur au moins une différence entre une partie imaginaire de valeur propre (Im(Àj)) évaluée dans le domaine temporel et une pulsation fondamentale (co). 7. Procédé selon l'une des revendications précédentes, caractérisé en ce qu'il comprend une étape (106) d'émission itérative ordonnancée pour réitérer lesétapes (101) à (104) en modifiant le mode vibratoire courant de façon à minimiser la fonction (f) tant que ladite norme est supérieure à la valeur de seuil (s). 8. Procédé selon l'une des revendications précédentes, caractérisé en ce qu'il comprend une étape (100) initiale ordonnancée pour évaluer des valeurs propres correspondant à des instabilités du système physique en équilibre statique, pour définir un mode vibratoire initial comprenant un vecteur initial de pulsations fondamentales (coo,j) en fonction des parties imaginaires des valeurs propres et un vecteur initial d'amplitudes (Ûo) des grandeurs physiques en fonction de vecteurs propres associés aux valeurs propres. 9. Procédé selon l'une des revendications précédentes, caractérisé en ce que le système physique est un système de frein pour véhicule. 10. Produit programme d'ordinateur chargeable directement dans la mémoire interne d'un calculateur numérique, comprenant des portions de code logiciel pour l'exécution des étapes de procédé selon l'une quelconque des revendications 1 à 7 lorsque ledit programme est exécuté sur un calculateur.