FR 2923045 A1 20090501 FR 0758728 A 20071031 Procédé de simulation par ordinateur d'un équilibre dynamique de système physique. L'invention concerne un procédé de simulation par ordinateur d'un équilibre dynamique de système physique. Le domaine de l'invention est celui du dimensionnement vibratoire d'un système physique, c'est-à-dire de l'évaluation et de l'ajustage de grandeurs physiques qui impactent le comportement du système en termes de phénomènes vibratoires. Dans la conception d'un système électrique, le dimensionnement consiste à déterminer les auto- inductances, les capacités et les résistances de différents composants d'un circuit. Dans la conception d'un système mécanique, le dimensionnement consiste à déterminer les masses, les raideurs, les coefficients d'amortissement et les formes de différents composants d'un assemblage. On connaît les méthodes de conception qui consistent a eviter ou amoindrir des phénomènes vibratoires gênants en modifiant la forme d'une pièce et en observant son comportement réel jusqu'à obtenir un comportement acceptable. Par exemple le document FR2415235A décrit une méthode de conception aboutissant à une forme de disque de frein qui réduit le crissement. La fabrication d'un prototype puis son essai en grandeur réelle, induit des coûts et des délais qui entravent une recherche tendue de performances sans cesse améliorées. L'état de la technique s'oriente depuis plusieurs années déjà vers des procédés de simulation assistés par ordinateur qui permettent d'améliorer le comportement, en particulier vibratoire, de systèmes physiques en évitant les contraintes de coût et de délais de fabrication et d'essais de prototypes matériels à différents stades d'évolution. Dans ce cadre, l'invention concerne un procédé de simulation qui utilise la méthode dite de la balance harmonique (HBM pour Harmonic Balance Method en anglais). La méthode de la balance harmonique est une méthode fréquentielle qu'il est déjà connu d'adapter de façon à traiter des problèmes non linéaires et d'obtenir rapidement une solution périodique approchée d'un système dynamique soumis à une excitation périodique ou pseudopériodique. Il est connu d'utiliser la méthode de la balance harmonique dans des systèmes stimulés par une vibration extérieure au système. Ainsi, le document W099/36872A qui divulgue une utilisation de la méthode de la balance harmonique pour simuler des circuits électroniques comportant une quantité importante d'éléments non linéaires, établit une distinction entre un signal d'entrée unique à une seule fréquence et plusieurs signaux d'entrée à différentes fréquences. Cependant, il existe des systèmes auto excités susceptibles de vibrer sans être stimulés par un signal périodique. C'est par exemple le cas d'un système de frein d'automobile qui génère des crissements lorsqu'il est actionné. Il existe de nombreux systèmes auto excités parmi lesquels certains génèrent des oscillations bénéfiques, comme c'est le cas des horloges, et d'autres génèrent des nuisances sonores voire destructives, comme il est préférable de ne pas les observer sur une structure de pont qui enjambe une rivière. Un besoin se fait donc sentir pour ne pas se limiter avec la méthode de balance harmonique au traitement des problèmes non linéaires dans les systèmes soumis à une excitation périodique ou pseudopériodique extérieure. L'invention a pour but de résoudre les problèmes vibratoires dans les systèmes auto excités. Un objet de l'invention est un procédé de simulation par ordinateur d'un équilibre dynamique de système physique soumis à des efforts non linéaires en générant un vecteur multi harmoniques qui minimise une fonction dynamique du système dans un domaine fréquentiel. De façon remarquable, le procédé comprend une étape d'évaluation dans laquelle un critère de minimisation porte sur au moins un paramètre associé à ladite fonction dynamique dans le domaine temporel. Cette étape d'évaluation permet non seulement l'estimation de la solution non-linéaire stationnaire du système sujet à une instabilité de flottement, mais aussi la détermination de la fréquence fondamentale de cette solution non-linéaire. Particulièrement, le paramètre est une valeur propre complexe de matrice de gain appliquée à un vecteur d'état associé dans le domaine temporel au vecteur multi harmonique. Avantageusement, le procédé comprend une étape de test dans laquelle une convergence vers le vecteur multi harmonique qui minimise la fonction dynamique, est détectée lorsque la valeur propre complexe a une partie réelle nulle. Avantageusement aussi, le procédé comprend une étape de descente vers le vecteur multi harmonique qui minimise la fonction dynamique dans laquelle une direction de descente est déterminée à partir des dérivées de la fonction dynamique portant sur au moins le paramètre associé à la fonction dynamique dans le domaine temporel. Plus particulièrement, le procédé comprend une étape de calcul de matrice de gain dans laquelle la matrice de gain est obtenue en additionnant une matrice d'état dynamique et une matrice jacobienne des efforts non linéaires massiques dans le domaine temporel et/ou une étape de calcul des efforts non linéaires dans laquelle une expression des efforts non linéaires est générée dans le domaine temporel puis convertie dans le domaine fréquentiel. L'invention sera mieux comprise, et d'autres buts, caractéristiques, détails et avantages de celle-ci apparaîtront plus clairement dans la description explicative qui va suivre faite en 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 montre des étapes du procédé selon l'invention pour simuler des comportements vibratoires stimulés ; la figure 2 montre des étapes de procédé selon l'invention pour simuler des comportements vibratoires 15 auto excités. En référence à la figure 1, on décrit une mise en oeuvre de la méthode de la balance harmonique appliquée à un problème mécanique numérique écrit selon la méthode des éléments finis. De façon connue, un maillage de 20 pièces mécaniques du système à simuler, définit un vecteur de grandeurs physiques U dont chaque composante qualifie un noeud du maillage. Pour simuler le comportement vibratoire du système, chaque composante du vecteur U représente un déplacement de noeud. Une 25 composante peut à son tour elle-même constituer un vecteur, par exemple, à trois sous composantes dans un espace tridimensionnel. Une dérivée par rapport au temps U du vecteur U représente une vitesse de chaque noeud correspondant et une dérivée par rapport au temps U du 30 vecteur U représente une accélération de chaque noeud correspondant. Le comportement vibratoire du système mécanique est modélisé par l'équation de mouvement exprimée sous la forme matricielle suivante : 35 MU + CU + KU + Fni = Fext (1) Dans laquelle : M est une matrice de masses C est une matrice d'amortissement K est une matrice de raideur Fext est un vecteur des forces extérieures Fn1 est un vecteur qui représente les forces non linéaires Le membre de gauche de l'équation (1) comporte les termes linéaires de la modélisation par éléments finis MU+CU+KU et l'effort distinction est nécessaire non linéaire Fn1. Cette car les termes linéaires et non linéaires ne sont pas traités de la même façon. La méthode de la balance harmonique repose sur une hypothèse physique selon laquelle, à une excitation périodique ou pseudo périodique, tout système va répondre de manière périodique. Ainsi, la solution se présente sous la forme d'un développement en série de Fourier : k=+oo k=+oo U(t) = k Uk cos(kcot) + k Usk sin(kcot) (2) k=0 k=1 Dans le cadre de l'invention la fréquence co constitue une inconnue du problème. L'accélération et la vitesse de l'équation (1) s'obtiennent en dérivant la série de Fourier par rapport au temps. Selon la méthode de la balance harmonique, on ne conserve que les premières harmoniques de (2) jusqu'à un rang N prédéterminé. Ainsi, pour un système physique dont l'essentiel de la dynamique est régi par ces harmoniques, la solution approchée peut s'écrire : k=Nh k=Nh U(t) Uapproché (t) = Uk cos(kcwt)+ U: sin(kcwt) ( 3 ) k=0 k=1 On rappelle que dans le cas d'un système purement 35 linéaire, on ne retient que la première harmonique. En effet, pour une excitation à une pulsation i la réponse a la même pulsation. Ce sont les efforts non linéaires Fn1 qui introduisent des couplages entre les différents coefficients de la série de Fourier. En remplaçant la solution réelle U(t) par la solution approchée Uapproché (t) dans l'équation du mouvement (1) on obtient un résidu RNh(t) dont l'expression est la suivante : RNh (t) = [(K ù k2cn 2M)U; + (kcwC)Uk ]cos(kcnt) + [(K ù k2cn 2M)Uk ù (kcwC)U; lsin(kcnt) k=O k=1 + FNL (t, U, U) - Fext (t) (4) En projetant le résidu sur la base des sinus et cosinus qui est une base orthonormée de l'espace des 15 fonctions continues, et en posant le vecteur multi harmoniques U tel que : LUcT CT ST T s T U = O U1 U1 UNä UNn 20 L'équation du mouvement (1) devient : AU + FNL = Fext Avec K 0 0 0 0 0 0 0 0 0 0 0 0 Ah (5) A= 25 (6) (7) [ù (kw)2M + K kwC ùkwC ù(kw)2M+K (8) où Ahk est la matrice de raideur dynamique associée à la ]dème harmonique. L'équation de mouvement (6) s'exprime dans le domaine fréquentiel. Cette dernière ne comprend que les coefficients de Fourier des degrés de liberté du système. Le principe de la balance harmonique est de résoudre ce système d'équations en équilibrant terme à terme les coefficients de Fourier de la solution approchée. Contrairement à la résolution des systèmes linéaires qui ne nécessitent pas de dispositions particulières, la résolution d'un système d'équations non linéaires requiert des méthodes mathématiques particulières. Dans le cas présent, cela revient à minimiser la fonction coût suivante : Î(U) = AU + FNL ù F'exc =0 (9) 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é, on rappelle simplement les éléments suivants. De nombreuses méthodes d'optimisation locale dérivent de la méthode de Newton dans laquelle, à partir d'un point initial sur la surface définie par les équations on se déplace de manière logique vers un minimum local de la fonction, c'est-à-dire en suivant la Ah,k = pour k =1...Nh pente de la fonction, comme un skieur qui descendrait la piste vers le télésiège. Tout comme le skieur qui descend les pistes en ligne droite 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 (Matlab, Scilab) ainsi que dans les librairies mathématiques gratuites sur Internet (GNU Scientific Library, Minpack, Levmar). La détermination des coefficients de Fourier des forces non linéaires est une tâche qui peut s'avérer difficile. Il est possible de déterminer directement les coefficients de Fourier des efforts non linéaires. Ces coefficients dépendent évidemment de U;+1 car ils sont non linéaires. Toutefois le travail analytique que demande cette technique est insurmontable et ne peut s'utiliser que pour des systèmes avec quelques degrés de liberté et un nombre d'harmoniques Ni, très faible (en général Si l'indice i est strictement positif, il provient d'une étape 106 dans laquelle il a été incrémenté et dans laquelle le vecteur U; courant est pris égal au vecteur U; qui avait été utilisé dans une exécution précédente de l'étape 104. L'étape 106 est exécutée chaque fois que l'étape 105 conclut négativement au test de convergence. Si l'indice i est nul, il provient de l'étape 100 dans laquelle il a été initialisé à zéro et dans laquelle le vecteur U; courant est initialisé à un vecteur U0. Dans l'étape 102, l'ordinateur calcule le vecteur 35 U;+1 au moyen de la formule : U,+1 := U, + a, s; La direction de descente s, est calculée dans une étape 101. Dans l'étape 101, l'ordinateur calcule la direction de descente s, à partir des dérivées de la fonction f(U,). Si i est strictement positif, la fonction f(U,) est celle déterminée dans une exécution précédente de l'étape 104. Si i est nul, la fonction f(U,) est une fonction f(Uo) initialisée dans l'étape 100. Le procédé présenté précédemment est bien adapté dans le cas des systèmes soumis à une excitation périodique. La figure 2 illustre les adaptations apportées au procédé lorsque le système est auto excité. Un système autonome est un système dynamique d'évolution temporelle continue qui a une indépendance explicite du temps t : X(t)=F(X) (10) Le système autonome ci-dessus présente plusieurs types de solutions dont deux principales qui nous intéressent, les points fixes et les cycles limites. En ce qui concerne les points fixes, on reprend l'équation (10) appliquée à un système autonome non linéaire : X(t)=F(X) (11) Un état d'équilibre est caractérisé par la relation suivante : F(Xe)=O (12) Toute solution Xe vérifiant cette relation est appelée position d'équilibre, point singulier, point fixe 35 ou encore solution stationnaire. Dans le cadre, par exemple, d'un système de freinage de véhicule automobile, la solution correspond30 au chargement statique des plaquettes de frein sur le disque dû à l'action du conducteur sur la pédale de frein. En ce qui concerne les cycles limites, on considère le système autonome non linéaire vu précédemment. Une solution périodique X(t) est une solution telle qu'il existe un réel T pour lequel : Quelque soit t, X(t +T) = X(t) et X(t +T) ~ X(t) pour 0 T est la période de la solution. La représentation d'une telle solution dans le plan de phase donne une trajectoire fermée appelée cycle limite. Les cycles limites représentent les états stationnaires d'oscillations. Le crissement de frein est un cycle limite, c'est-à-dire que c'est la solution périodique stationnaire du système de frein. Dans le cas de systèmes non linéaires tels que le système de freinage, la théorie des cycles limites est en relation avec la bifurcation de Hopf, ce qui est expliqué plus loin dans la description. Il y a différents sens que l'on peut donner quant à la stabilité d'une solution. Pour des raisons de simplifications, seules deux définitions de stabilité concernant le système sont étudiées. La première définition est celle au sens de Lyapunov qui s'énonce de la façon suivante : La solution X(t) est stable si pour tout e ,>O, il existe ô (E3)>O tel que toute solution Y(t) vérifiant Y(to)-X(to) 13 La solution X(t) est asymptotiquement stable si elle est stable au sens de Lyapunov et si en plus on a lim Y(t) ù X(t) = 0 t-+oo Revenant sur les points fixes, on considère l'équation de système autonome non linéaire (10). Dans un cadre de linéarisation du système de départ, considérant un point singulier Xe du système non linéaire vérifiant l'équation (12), une perturbation Y génère une valeur X telle que : X(t) = Xe (t) + Y(t) (15) En substituant (15) dans (10) on obtient : 15 Y=F(Xe+Y) (16) La position d'équilibre X=Xe est maintenant devenue Y=0. En supposant que la fonction F est de classe Cl, son développement en série de Taylor en Xe à l'ordre 1 donne : Y = F(Xe)+DXF(XQ)Y+o(Y2) (17) Et donc Y DxF(XQ)Y = JFY (18) 25 La recherche d'une solution de l'équation (18) équivaut sous forme linéarisée à la recherche d'une solution de l'équation (10) au niveau du point fixe X=Xe. L'examen des valeurs propres d'une matrice jacobienne JF 30 permet de conclure quant à la stabilité du système autour du point fixe. On rappelle que la matrice jacobienne d'une fonction vectorielle F d'un vecteur X, F(X)= ~f1(xir...,Xk,...,xn),...,f (x1,...,Xk,...,xn),...,fn(x1,...,Xk,...,Xn) }T en j=n k=n X=Xe est la matrice des dérivées partielles JF = af' (Xe) axk j=1 k=1 20 Soit X (j=1, n) les valeurs propres de la matrice JF. Considérant la partie réelle de chaque valeur propre, d'après le théorème de Lyapunov, il existe trois cas possibles. Dans un premier cas, Vj,Re(X1) Dans un second cas, il existe au moins un indice i tel que Re(X1) > 0 alors Y=0 est un point fixe asymptotiquement instable de Y = JFY. Cela signifie qu'une perturbation du point fixe entraîne des oscillations bornées du système. Dans un troisième cas, si Re(X1) 0. Le système est asymptotiquement instable. Dans la deuxième zone, une perturbation du point fixe entraîne des oscillations de période T =cw/2n qui tendent vers un cycle limite lorsque t-* +oo. Ainsi, la bifurcation de Hopf correspond au passage d'un état d'équilibre statique à un état d'équilibre dynamique. Après un point de bifurcation de Hopf, deux solutions peuvent coexister : une solution d'équilibre statique et une solution d'équilibre dynamique. Le passage de la solution d'équilibre statique à la solution d'équilibre dynamique requiert une modification des paramètres du système tel ü ou une perturbation des conditions initiales. Lorsqu'une intégration temporelle est réalisée et que le système se situe après le point de bifurcation, la non perturbation des conditions initiales donne une solution statique qui reproduit les conditions initiales. Par contre une perturbation sur un paramètre ou bien sur les conditions initiales provoque une solution périodique stationnaire. Le cas de figure qui permet d'atteindre le cycle limite, avec une intégration temporelle, ne fonctionne pas avec l'utilisation de la HBM telle qu'elle est mise en oeuvre sans utiliser l'invention. Pour rappel, le principe de la HBM est de trouver un minimum de l'équation dynamique dans le domaine fréquentiel. Après le point de bifurcation, on a vu qu'il coexiste deux solutions possibles. Ainsi, la méthode HBM a le choix entre deux minima, l'un correspondant à l'équilibre statique (trivial qui correspond à une solution exacte de l'équation (9), mais définissant un équilibre instable) et l'autre correspondant à l'équilibre dynamique (cycle limite définissant la solution stationnaire stable plus particulièrement intéressant dans le cadre de l'invention). A noter que la solution statique qui correspond à la solution exacte du problème à minimiser sera toujours la solution trouvée par la méthode HBM classique quelque soit les conditions initiales, cette solution statique correspondant à celle où tous les coefficients de la série de Fourier sont à 0 sauf ceux décrivant les termes constants. Ainsi, il est nécessaire à ce stade de contraindre l'équation (9) afin de forcer un mécanisme solveur à aller sur la solution périodique dont la pulsation fondamentale i est une inconnue du problème. Le procédé qui fait l'objet de l'invention, comprend avantageusement des moyens pour forcer le mécanisme solveur à aller sur la solution périodique en recherchant la fréquence fondamentale associée. Ces moyens sont décrits à présent dans le cadre d'un exemple de mise en oeuvre sur un système de freinage de véhicule automobile. En reprenant l'équation (1) et en l'écrivant dans un espace d'état, on obtient : IÙ..1 [A] U + 0 + 0(23) U U M Fext ûM-1F 1 Avec 0 I [A] est une matrice d'état du système dynamique appliquée au vecteur d'état X qui décrit le système dans l'espace des phases en posant : U X= iJ Soit X = AX + M 'FExT -M 'FNL (X) Si on applique la méthode de perturbation vue plus haut au système ci-dessus, on obtient le système 25 linéarisé suivant : @ AY+Dx[-M-lFn1(Xp)]Y=AY+Jä1Y (27) Avec Jn1 la matrice jacobienne de l'effort non 30 linéaire multiplié par l'inverse de la matrice des masses. L'étude des valeurs propres de (A+Jnl) permet de conclure quant à l'état de la solution calculée Xp(t) lors de l'optimisation. Parmi les valeurs propres [Al= M- 'K M 1C û (24) 20 (25) (26) Xj (j=1, n) de la matrice (A+Jw), il existe h tel que ?L'h =instable , c'est-à-dire qu'une des valeurs propres correspond au mode instable du système qui génère les oscillations. Ainsi, si Re(Xinstable)=0 alors la solution calculée Xp(t) est la solution du système et donc le système a atteint l'équilibre dynamique, à savoir une vibration à fréquence et amplitude constante. Néanmoins, si la valeur de la partie réelle est différente de 0, alors la variable temporelle Xp(t) n'est pas solution du système dynamique. Dans ce cas, Xp(t) correspond à un état transitoire. Dans ce problème d'optimisation où deux solutions coexistent, l'intérêt de rajouter une équation supplémentaire est de contraindre le mécanisme solveur à atteindre la solution périodique recherchée. Ainsi, la fonction à minimiser du problème d'optimisation (9) s'écrit dorénavant sous la forme f(U,w) = f(Ù,w) o f2(U,w) o avec fl(U,O))=A(0))U+FNL(U)ùFext (28) et f2 (U,(0) = Re(Xinstable ) C'est le rajout de cette équation supplémentaire qui fait l'originalité du procédé selon l'invention et permet la détermination non seulement de la solution stationnaire non-linéaire, mais aussi de la fréquence fondamentale associée. Ainsi, lorsque l'optimisation sera terminée et si le mécanisme solveur a convergé, la solution trouvée correspondra au cycle limite du système dynamique. Le procédé décrit à présent, en référence à la figure 2, comprend une étape finale 110 qui donne une valeur de la pulsation inconnue i et un vecteur multi harmonique Usumtion à n•l• (2m+1) composantes où m désigne un nombre d'harmoniques d'une pulsation de base c.) et n désigne un nombre de noeuds d'une modélisation par éléments finis du système mécanique pour lequel les composantes du vecteur multi harmonique Usomtion sont une solution d'un équilibre dynamique du système mécanique dans les conditions physiques ou environnementales dans lesquelles le système mécanique est étudié. Dans un espace physique à trois dimensions tel qu'on le considère habituellement en mécanique newtonienne, le nombre 1 de degrés de liberté de chaque noeud est égal à 3. Les n•l composantes du vecteur Usomtion associées à l'harmonique d'ordre zéro, correspondent à une partie statique de l'équilibre dynamique du système mécanique. Les autres composantes du vecteur Usomtion associées aux harmoniques d'ordre 1 et supérieures, correspondent à une partie vibratoire de l'équilibre dynamique du système mécanique. Le vecteur Usomtion est égal à un vecteur U;+1 détecté dans une étape de test 115. Dans l'étape 115, un ordinateur teste si la valeur du vecteur Ui+1 obtenu après i+1 exécutions d'une étape d'évaluation 109, a convergé vers une solution recherchée. La convergence du vecteur U;+1 est considérée atteinte dès que l'ordinateur détecte une composante nulle dans un vecteur f2(U,+1,co) calculé dans l'étape 109. Dans l'étape 109, l'ordinateur calcule un vecteur f (U,+1,co) qui comprend un vecteur f (U,+l,co) et le vecteur f2(Ui+1,c)) . Le vecteur f (Ui+1,co) est calculé en multipliant le vecteur U;+1 par une matrice de raideur dynamique A(w) et en combinant le résultat obtenu à deux vecteurs d'efforts FjUi+1) et Fe, dans le domaine fréquentiel. Un vecteur d'efforts non linéaires Fni(U,+1) est ajouté au résultat obtenu et un vecteur d'efforts extérieurs Fe.t est retranché au résultat obtenu. Le vecteur d'efforts non linéaires Fni(U,+1) qui dépend du vecteur multi harmoniques de déplacements U;+1, est préalablement calculé dans une étape 113 pour chaque exécution de l'étape 109. Le vecteur d'efforts extérieurs Fe.t qui ne dépend pas du vecteur multi harmoniques de déplacements U;+1, est déterminé une fois pour toutes dans une étape d'initialisation 110 décrite plus loin. Par contre la matrice de raideur dynamique A(w) qui dépend de la fréquence fondamentale ici inconnue doit être déterminée à chaque évaluation de lors de f2(Ui+1,{o) l'optimisation. Le vecteur comprend la partie réelle du mode instable calculé préalablement dans une étape 108 pour chaque exécution de l'étape 109. Dans l'étape 108, l'ordinateur calcule une matrice de gain J en additionnant une matrice d'état dynamique A et une matrice jacobienne Jn1 des efforts non linéaires massiques M-1FNL(X) dans le domaine temporel. La matrice d'état dynamique A et la matrice de masse M sont des données intrinsèques du système mécanique qui ne dépendent pas du vecteur multi harmoniques de déplacements U;+1. Par conséquent ces deux matrices sont déterminées une fois pour toutes dans l'étape d'initialisation 110. Les efforts non linéaires FNL(X) dépendent d'un vecteur d'état X qui comprend les vecteurs de déplacement U et U associés dans le domaine temporel au vecteur multi harmoniques de déplacements U;+1. Par conséquent, le vecteur des efforts non linéaires FNL(X) est préalablement calculé dans l'étape 113 pour chaque exécution de l'étape 108. Le calcul des valeurs propres À de la matrice de gain J est effectué par une méthode classique. Dans l'étape 113, l'ordinateur effectue une transformée inverse rapide FFT-1 du vecteur U;+1 de façon à obtenir les vecteurs de déplacement U et U associés dans le domaine temporel. Il est plus facile de traiter les non linéarités dans le domaine temporel. On obtient ainsi un vecteur déplacement Ui+l = FFT 1 (U;+1) qui permet de faire le calcul des efforts non-linéaires en faisant intervenir des raideurs non linéaires continues de type cubiques et une gestion de la perte de contact permettant de représenter l'évolution des contacts des pièces mécaniques telles que par exemple des plaquettes de frein sur un disque soumises à une pression de commande de freinage. On obtient alors dans le domaine temporel un vecteur de forces non linéaires Fn1 (Ui+1) sous la forme : Fn1 (Ui+1) - 8 contactkUi+13 avec 6contact égal à 1 si il y a contact et 0 dans le cas de non contact. Une transformée rapide FFT est appliquée au vecteur Fn1 (Ui+1) de façon à obtenir les efforts non linéaires Fn1(U,+1) dans le domaine fréquentiel sous la forme : Fnl(Ui+1) = FFT (Fn1 (Ui+1) Le vecteur Ui+1 utilisé dans les étapes 113, 109 et 115, est préalablement calculé pour l'exécution de chacune de ces étapes dans une étape 112 par une méthode de descente à partir d'un vecteur U; précédemment connu. Si l'indice i est strictement positif, il provient d'une étape 116 dans laquelle il a été incrémenté et dans laquelle le vecteur U; courant est pris égal au vecteur U; qui avait été utilisé dans une exécution précédente de l'étape 109. L'étape 116 est exécutée chaque fois que l'étape 115 conclue négativement au test de convergence. 30 Si l'indice i est laquelle il a été vecteur U; courant Dans l'étape Ui+1 au moyen de la Ui+1 =U,+a,s, La directionnul, il provient de l'étape 110 dans initialisé à zéro et dans laquelle le est initialisé à un vecteur Uo. 112, l'ordinateur calcule le vecteur formule . de descente s, est calculée dans une étape 111. Le coefficient multiplicateur ai peut être pris égal à une constante ao. Dans ce cas la constante ao est judicieusement choisie pas trop petite pour ne pas 35 progresser trop lentement vers la solution et pas trop grande pour ne pas dépasser la solution alternativement dans un sens et dans un autre. Le coefficient multiplicateur ai peut aussi être multiplié par un coefficient p inférieur à 1 à chaque itération de façon à diminuer progressivement au fur et à mesure qu'on se rapproche de la solution. Dans l'étape 111, l'ordinateur calcule la direction de descente s, à partir des dérivées de la fonction f(U,,co) . Si i est strictement positif, la fonction f(U,,co) est celle déterminée dans une exécution précédente de l'étape 109. Si i est nul, la fonction f(U,,co) est une fonction f(U0,c0o) initialisée dans l'étape 110. On note que la direction de descente s, combine le vecteur f (U,,co) dans le domaine fréquentiel avec le vecteur f2(Ui+14o) dans le domaine temporel. L'étape d'initialisation 110 prend en compte une modélisation par éléments finis d'un système mécanique. En particulier, le procédé rapatrie de la modélisation par éléments finis, les matrices de masse M, d'amortissement C et de raideur K de façon à calculer la matrice A qui est une constante au regard du traitement des instabilités par le procédé : 0 I ùM-1K ùM-1C Un déplacement { ux, uy, uZ } est associé à chaque noeud du modèle par éléments finis. L'ensemble des déplacements associés est regroupé dans un vecteur de déplacement U à 3n composantes où n est le nombre de noeuds du modèle. Au vecteur de déplacement U, on fait correspondre le vecteur multi harmonique U; où i est un indice d'itération qui sera utilisé pour exécuter les étapes suivantes jusqu'à la convergence vers une solution. Dans l'étape d'initialisation 110, l'indice i d'itération est initialisé à zéro. Le vecteur multi harmonique U; a 3n(2m+1) composantes où m désigne le nombre d'harmoniques d'une pulsation de base w. Il est initialisé à une valeur initiale U0. Le procédé converge avec, par exemple, toutes les composantes nulles dans le vecteur initiale U0. Cependant, la convergence risque [A] _ d'être lente. Pour améliorer la rapidité de convergence du procédé, on prend une valeur initiale égale au déterminant de la matrice A et on donne une valeur initiale pour le vecteur Uo qui définit la solution divergente dans le domaine temporel, solution approchée par le mode instable de la solution, et donc donné sous la forme Y(t) = p (Ive' + jet où définit le conjugué complexe de i avec Y(t) solution divergente du système i vecteur propre du mode instable du système linéarisé X valeur propre complexe du mode instable du système linéarisé p constante d'initialisation arbitrairement choisie 15 Dans l'étape 110, le vecteur d'efforts extérieurs Fe.t est déterminé à partir des conditions aux frontières du modèle telles que, par exemple, dans le cas d'un système de freinage automobile, les forces d'appui sur 20 les plaquettes en contact avec le disque et le couple exercé sur le disque par le moyeu des roues. Dans le cadre d'une simulation de système autonome, le vecteur Fe.t est limité à sa composante d'harmonique zéro si on considère les vibrations engendrées par une application 25 de forces statiques. Dans l'étape 110, la matrice de raideur dynamique A est initialisée à la valeur : K 0 0 0 A= 0 Am 0 0 0 0 0 0 0 0 Ah Avec [ù (kw)2M + K kwC ùkwC ù(kw)2M+K Ah,k pour k =1...Nh 30 Dans l'étape 110, on définit aussi une fonction f(U,,co) initialisée à f(Uo,coo) en prenant pour ses composantes dans le domaine fréquentiel le vecteur : f (Uo,{oo) = A(ao)Uo ùFexr et en prenant pour ses composantes dans le domaine temporel le vecteur : f2(Uo,cwo)=Re(a,1) où Xj est la valeur propre de la matrice A correspondant au mode instable. L'enseignement fait ici dans le cadre d'un système mécanique est aisément transposable à d'autres types de systèmes physiques. Par exemple pour un système électrique, il suffit de remplacer la matrice des masses par une matrice des inductances, la matrice des amortissements par une matrice des résistances et la matrice des raideurs par une matrice des capacitances. Dans un système électrique, le vecteur U désigne une charge électrique et sa dérivée par rapport au temps un courant électrique. REVENDICATIONS 1. Procédé de simulation par ordinateur d'un équilibre dynamique de système physique soumis à des efforts non linéaires en générant un vecteur multi harmonique qui minimise une fonction dynamique du système à fréquence inconnue dans un domaine fréquentiel, comprenant : - une étape d'évaluation (109) dans laquelle un critère de minimisation porte sur au moins un paramètre associé à ladite fonction dynamique dans le domaine temporel et permet d'estimer la solution dynamique en terme de fréquence et d'amplitudes vibratoires 2. Procédé selon la revendication 1, dans lequel ledit paramètre est une valeur propre complexe de matrice de gain dont on cherche à minimiser la partie réelle en minimisant ladite fonction dynamique. 3. Procédé selon la revendication 2, comprenant : une étape de test (115) dans laquelle une convergence est détectée lorsque le vecteur mufti• harmonique minimise la fonction dynamique et lorsque la valeur propre complexe a une partie réelle nulle. 4. Procédé selon l'une des revendications précédentes, comprenant : - une étape de descente (111) vers le vecteur multi harmonique qui minimise la fonction dynamique dans laquelle une direction de descente est déterminée à partir des dérivées de la fonction dynamique portant sur au moins ledit paramètre associé à ladite fonction dynamique dans le domaine temporel. 5. Procédé selon l'une des revendications 2 ou 3, comprenant . 26 - une étape de calcul de matrice de gain (108) dans laquelle la matrice de gain est obtenue en additionnant une matrice d'état dynamique et une matrice jacobienne des efforts non linéaires massiques dans le domaine temporel. 6. Procédé selon l'une des revendications précédentes, comprenant : - une étape de calcul des efforts non linéaires (113) dans laquelle une expression des efforts non linéaires est générée dans le domaine temporel puis convertie dans le domaine fréquentiel. 7. Procédé selon l'une des revendications précédentes, dans lequel le système physique est un système de frein, notamment pour véhicule automobile, système auto-entretenu par frottement, la fonction dynamique du système s'exprimant par la formule suivante : 1+{ 01 (23) où U U M F,', _m 25 [A] est une matrice d'état du système appliquée à un vecteur d'état X qui décrit le système dans un espace des phases M est une matrice de masses F~ 8. Procédé selon la revendication 7, dans lequel le 30 vecteur d'efforts non linéaires fait intervenir des raideurs non linéaires continues de type polynomiales et une gestion de la perte de contact, de nature discontinue, permettant de représenter l'évolution des contacts au niveau d'un disque et de plaquettes de frein. 35 27 9. Procédé selon la revendication 7 ou 8, dans lequel le vecteur d'efforts extérieurs est représentatif de forces d'appui sur des plaquettes de frein en contact avec un disque et d'un couple exercé sur le disque par un moyeu de roues.