FR 2966950 A1 20120504 FR 1058878 A 20101028 PROCEDE DE SIMULATION DE COMPORTEMENT DE SYSTEMES PHYSIQUES, EN PARTICULIER DE CIRCUITS HYDROMECANIQUES INTRODUCTION L'invention concerne un procédé de simulation comportemental de systèmes physiques en temps réel, en particulier de circuits hydromécaniques complexes, tout en préservant la structure topologique - respectivement l'architecture hydromécanique - de ces circuits. Une application vise la commande hydromécanique d'une colonne de direction d'un aéronef. DOMAINE TECHNIQUE La présente invention se rapporte au domaine de la prévision du comportement des systèmes physiques, en particulier des systèmes hydromécaniques complexes, couvrant un spectre de basses fréquences de 0 à 20 Hz, afin de fournir une information exploitable en temps réel sur le comportement des paramètres dynamiques fondamentaux, par exemple la pression et le débit des fluides dans un circuit hydromécanique. Il s'agit alors de modéliser les paramètres dynamiques fondamentaux, de type « flux » - débit, couple, force, intensité électrique - et « potentiel » - pression, vitesse de rotation, vitesse linéaire, tension électrique -, selon le domaine technique afin de simuler en temps réel la mise en condition des fluides avec un niveau de performance élevé c'est-à-dire, dans le domaine de la simulation, une représentativité réaliste, notamment en cas de défaillance du circuit. ETAT DE LA TECHNIQUE ANTERIEURE ET PROBLEMATIQUE Dans le but de simuler en temps réel le comportement de certains systèmes physiques complexes, les procédés de simulation actuels remplacent l'ensemble du système modélisé par un filtre du premier ordre. Ces solutions, développées par exemple par le logiciel de simulation AMESIM de la société LMS, ne permettent pas de préserver la représentation du système en temps réel : soit les solutions obtenues ne sont pas compatibles avec une exploitation en temps réel, soit les simplifications opérées ne respectent pas la topologie du système. En effet, une gestion du comportement en temps réel des modèles hydromécaniques, par exemple par simulation d'architectures d'aéronefs en vol, peut être définie par une connaissance des paramètres fondamentaux dans un intervalle de temps de 10 ms par étape de simulation, soit jusqu'à 100Hz dans le domaine fréquentiel. Concrètement, pour conserver une marge suffisante permettant d'éviter les problèmes d'instabilité numérique, l'intervalle de fréquences des phénomènes physiques qu'il convient de représenter se situe en fait dans un intervalle de fréquences plus restreint, inférieures ou égales à 20 Hz, en divisant par 5 l'intervalle théorique de 100Hz. Cependant, une durée de 10 ms par étape sera retenue tout en conservant cette marge lorsque la simulation fournit des résultats encore satisfaisants - au regard des paramètres physiques en jeu - jusqu'à 200 ms par étape de simulation. Dans ces conditions, les phénomènes de basses fréquences, représentatifs des phénomènes physiques fondamentaux, restent pris en compte alors que les phénomènes de vibration ou de dynamique à hautes fréquences sont filtrés. Ainsi, d'une part, les simulations connues sont incapables de filtrer les phénomènes hautes fréquences et ne permettent donc pas une modélisation représentative du comportement : elles restent liées au principe de causalité intrinsèque qui oblige à résoudre implicitement, dans tous les cas, des équations dynamiques dans le domaine des hautes fréquences, à savoir des équations de compressibilité en hydraulique et des équations de dynamique en mécanique, aux noeuds du circuit. Ainsi, la dynamique des moyens de régulation, tels que des servo-valves, est généralement modélisée par un filtre du second ordre. Dans ces conditions, de telles gestions tiennent compte de la raideur inertielle de l'architecture du circuit, qui génère des hautes fréquences, en général supérieures à 100 Hz. D'autre part, ces simulations ne permettent pas d'atteindre une gestion en temps réel car le comportement hydraulique reste difficile à représenter du fait de la non-linéarité de la relation de base entre les paramètres dynamiques : par exemple, la relation de base Q(P) entre le débit Q et la variation de pressions OP d'un fluide dans un restricteur de flux est fournie par l'expression Q = k(P)-/IOPI, k étant le gain qui dépend de la pression avec une forte non-linéarité pour les faibles pressions. Cette non-linéarité oblige à utiliser des algorithmes de résolution non linéaires dont la convergence n'est pas assurée dans un laps de temps réduit : la contrainte de gestion en temps réel n'est alors pas respectée. OBJET DE L'INVENTION L'invention vise à pallier ces insuffisances en permettant une simulation du comportement d'un système physique complexe en temps réel - avec par exemple une durée de simulation par étape de 10 ms - suffisamment représentative par une modélisation simplifiée des composants, tout en préservant la topologie du système et ainsi permettre une représentativité efficace des cas de défaillance. Pour ce faire, l'invention se base sur une approche nodale à partir de liens graphiques - et non de liens de causalité - entre composants du circuit pour déterminer les relations entre paramètres aux noeuds, puis procéder à une linéarisation par une approche par corde et à une résolution des équations linéaires. Plus précisément, la présente invention a pour objet un procédé de simulation de comportement d'un système physique complexe consistant à créer une librairie de composants, à assembler des composants de cette librairie comportant au moins un composant de comportement non linéaire pour modéliser une architecture de circuit par interconnexion des ports d'entrée des composants à travers des noeuds d'accès, chaque port étant connecté à un seul noeud, à caractériser chaque port de chaque composant par un couple de paramètres potentiel/flux, à linéariser le comportement du (des) composant(s) non linéaire(s) selon une approximation par corde et à établir une relation de détermination du flux à partir de la variation de potentiel à chaque port des composants de comportements linéaires et linéarisé(s), le potentiel de chaque port étant assimilé à celui de son noeud de connexion, à constituer un système d'équations linéaires et linéarisée(s) entre les paramètres en appliquant la loi des noeuds en chaque noeud et à résoudre ce système d'équations de manière itérative pour déterminer les valeurs des paramètres à chaque noeud. Dans le cas de circuits hydromécaniques : - le circuit comporte au moins un composant non linéaire intemporel, c'est-à-dire dont le comportement n'est pas lié par une variation des paramètres en fonction du temps, en particulier par une dérivée de paramètre en fonction du temps ; - les paramètres de potentiel sont notamment : la pression, la vitesse angulaire ou la vitesse linéaire, et les paramètres de flux : le débit, le couple ou la force; - les composants non linéaires intemporels plus particulièrement concernés, mais non exclusivement, sont des générateurs 20 de pression, des servovalves, des valves (de sécurité, de contrôle, anti- shimmy, etc.), des restricteurs de flux et des vérins ; - afin d'améliorer les performances du procédé en représentativité, sans nuire à la compatibilité au temps réel, la loi des noeuds peut avantageusement faire intervenir un facteur de compressibilité temporel 25 à linéariser selon une approche par linéarisation tangentielle avec schéma implicite ; cette approche est formalisée par un schéma implicite substituant la dérivée partielle de potentiel en fonction du temps par une variation de ce potentiel par rapport à la variation de temps finie correspondante ; cette prise en compte de la compressibilité des fluides permet d'intégrer un facteur de 30 compressibilité réaliste et ainsi de rendre compte d'adaptations transitoires qui se traduiraient sinon en simulation par des situations de blocage. 15 Selon des modes de mise en oeuvre plus particuliers : - le système d'équations linéaires et linéarisée(s) entre le paramètre de flux et le paramètre de potentiel pour chaque noeud autour d'un composant Cj, à partir des variations du paramètre de potentiel à un instant antérieur selon l'approximation par corde pour le(s) comportement(s) non linéaire(s), est transposé en mode matriciel par formation d'une matrice de modélisation [Al et d'un vecteur complémentaire [BI, puis pour l'ensemble des composants Cj du circuit par intégration des matrices [Al dans une matrice de modélisation globale [A] et des vecteurs complémentaires [BI dans un vecteur global [B] selon une organisation définie par une matrice de topologie [N] du circuit ; - la topologie d'un circuit est modélisée en numérotant chaque composant, chaque port de chaque composant et chaque noeud, puis en positionnant les numérotations des noeuds en mode matriciel dans la matrice [N] de topologie en ligne selon la numérotation des composants et en colonne selon la numérotation des ports ; - un circuit étant modélisé, la simulation de son comportement par les valeurs de flux aux ports de chaque composant Cj à partir des valeurs de potentiel à chaque noeud par application de la loi des noeuds, est déterminée par des matrices [Ain et vecteurs [B]n de contribution des composants du circuit en chaque noeud « n » sont intégrées - par pavé - dans des matrices globales de circuit [A]c et [B]c, selon la topologie de position des noeuds propre au circuit, les valeurs de potentiels étant fournies par la résolution matricielle ([p] _ - [B]/[A])c, et la position des noeuds étant définie pour chaque composant en fonction des ports par la matrice de topologie [Mc, les valeurs de potentiel aux ports d'un composant «j » donné sont alors fournies par extraction des potentiels [pic correspondant aux noeuds identifiés dans la matrice de topologie [Mc, les valeurs de flux sont alors fournies par la relation flux potentiel. BREVE DESCRIPTION DES DESSINS D'autres caractéristiques et avantages de la présente invention apparaîtront à la lumière de la description qui suit concernant des exemples de réalisation non limitatifs, en référence aux dessins annexés qui représentent, respectivement : - en figures la et 1 b, un exemple d'application du procédé de l'invention à un restricteur de flux modèle et à un ensemble de trois restricteurs en série ; - en figure 2, un exemple d'application à un circuit élémentaire réduit à un seul composant, à savoir un générateur de pression ; et - en figure 3, un exemple d'application du procédé de l'invention sur un circuit de commande simplifié de direction hydromécanique de roulette d'avion. EXPOSE DETAILLE DE MODES DE REALISATION PARTICULIERS 15 La détermination des valeurs des pressions et des débits en fonction des pressions - plus généralement des paramètres de potentiels et de flux - est réalisée, dans un circuit hydromécanique comportant des composants non linéaires connectés entre eux par des conduits, en considérant la loi des noeuds à chaque noeud, et la linéarisation des relations 20 entre les paramètres de flux et de potentiel par une approximation par corde. Une résolution numérique dans une unité de calcul appropriée est effectuée par transcription matricielle, particulièrement avantageuse dans le cas d'un circuit complexe. Ces trois phases sont développées ci-après en partant de 25 portions de circuit comportant un seul composant de base : un restricteur de flux modèle, trois restricteurs en série et un générateur de pression. Le procédé est ensuite illustré sur un circuit hydromécanique de commande de direction hydromécanique, faisant intervenir des composants tels que des restricteurs et un générateur de pression, mais également d'autres 30 composants plus complexes. La figure 1 a représente schématiquement une portion élémentaire 1 de circuit entre deux noeuds Ni et Nj, comportant un restricteur 10 de flux 10, ayant deux ports d'entrée 12 et 14 de fluides, respectivement à des pressions P; et Pi, égales respectivement aux pressions PM et PNB aux noeuds Ni et Nj - selon le principe de l'invention qui veut que le potentiel de chaque port est assimilé à celui de son noeud de connexion - et selon des débits Q; et Qi. Un restricteur est un limitateur de flux, qui rend compte des variations de flux d'un conduit, telle que celles provoquées par une vanne ou un régulateur. Selon la loi des noeuds, la somme des débits est nulle (Q; + Qi = 0) dans le cas où les fluides sont considérés comme incompressibles et les débits ont alors pour valeur, comme rappelé en introduction : Qi = - Q; = K ~IoPI , OP étant égal à Pi - P; et K le gain du restricteur. Selon la convention usuelle des signes, les débits entrants sont positifs et les débits sortants négatifs. La linéarisation en utilisant l'approximation par corde passant par l'origine fournit alors les valeurs suivantes : Qi = -Qi = K OP V APant OPant désignant la variation de pressions OP à un instant antérieur. L'approximation consiste à substituer OPant à ./IoPI , ce qui revient à assimiler la courbe de variation Q = K-JIOPI à la corde qui la soutient telle que définie entre la valeur à l'origine de la courbe et la valeur de la courbe à l'instant antérieur. L'approximation par corde est plus performante que l'approximation tangentielle dans un calcul itératif à partir d'au moins deux relations flux/potentiel, tel que celui auquel la présente invention conduit. En effet, la substitution par la tangente fournit, à chaque étape de l'itération, des valeurs de convergence intermédiaires plus éloignées de la valeur d'équilibre recherchée que celles fournies par l'approximation par corde. En effet, l'approximation par corde sous-tend, en associant une valeur courante de chaque fonction à sa valeur à l'origine, les fonctions flux/potentiel qui présentent une concavité du type de celle la fonction racine carrée - fonction fondamentale dans les systèmes hydrauliques - alors que l'approximation tangentielle « redresse » la courbure de ces fonctions flux/potentiels. Ainsi, l'approximation par corde permet d'atteindre la valeur d'équilibre en un nombre de pas sensiblement inférieur à celui de l'approximation tangentielle, cette dernière pouvant même diverger du fait des conditions initiales. Des exemples de résolution en fonction de conditions initiales sont fournis ci-après. En référence à la figure 1 b, le premier exemple se rapporte à une portion de circuit 2 de circulation d'un fluide comportant un ensemble de trois restricteurs 10a, 10b et 10c en série, du même type que le restricteur selon la figure 1 a (respectivement de ports d'entrée 12a et 14a, 12b et 14b, 12c et 14c) chacun des restricteurs ayant le même gain intrinsèque K. Par exemple, une pression de charge PA de 10 bar (106 Pa) est fournie au port d'entrée 12a provenant d'un générateur 11 et une pression de retour PB à la bâche (réservoir 13) de 0 bar est recueillie au port de sortie 14a. Les pressions P1 et P2 aux noeuds N1 et N2 entre restricteurs sont déterminées à partir des relations issues de la loi des noeuds entre les restricteurs 10a et 10b d'une part, 10b et 10c d'autre part. Dans ces relations, les expressions des débits en fonction des pressions sont linéarisées, de manière similaire au traitement ci-dessus qui se rapportait à un restricteur unique. Le système d'équations aux noeuds N1 et N2 s'écrit alors respectivement : K1 (PA - P1) = K2 (P1 - P2) (I) K2 (P1 - P2) = K3 (P2 - PB) La linéarisation par l'approximation par corde introduit ainsi - dans les relations (I) - des gains de linéarisation K1, K2 et K3 qui dépendent du gain intrinsèque K et des valeurs de variations de pressions (PA - P1)ant, (P1 - P2)ant et (P2 - PB)ant à un instant antérieur, respectivement selon les expressions : K~ K K , K2 = = , K3 = (II) .\/(PA -P 1). ,I(Pl -P2 )ant ''/(P2 -PB ) ant Bien entendu, les valeurs des pressions PI et P2 à l'instant antérieur ne doivent pas annuler le dénominateur des expressions de K1, K2 et K3 ci-dessus. Si une valeur de PI ou P2 annulait l'un de ces dénominateurs, ce dernier serait remplacé par une valeur seuil sn, minimale prédéterminée de manière à ne pas bloquer les calculs. La résolution analytique du système I, avec - pour simplifier les expressions - une pression de retour PB nulle, fournit alors les solutions de pressions PI et P2 suivantes, en fonction des gains de linéarisation et de la pression de charge PA : KlK2 +KlK3 P~ PA KlK2 +KlK3+K2K3 P2= K1K2 PA KlK2 +KlK3+K2K3 En fonction des conditions initiales à l'instant to, se rapportant ici aux valeurs initiales de PA, PB, Pl et P2, des premières valeurs des gains de linéarisation sont calculées puis des premières valeurs des pressions PI et P2 aux noeuds. Ces premières valeurs sont réintroduites dans les expressions des gains pour former un calcul itératif entre les gains et les pressions. Les valeurs des gains et des pressions convergent vers des valeurs d'équilibre qui ne dépendent que du circuit 2, quelles que soient les conditions initiales, les valeurs initiales conditionnant le nombre de pas de calcul à effectuer avant d'atteindre ces valeurs d'équilibre. Ainsi, en partant des conditions initiales à l'instant to (PI) to = 2/3 PA et (P2) to = 1 /3 PA Le premier pas d'itération (correspondant à l'instant ti) fournit à partir des relations II et III : K~ = K2 = K3 = K P~ = 2/3 PA, P2 = 1 /3 PA Vl/3PA Un seul pas de calcul suffit dans ce cas particulier pour aboutir K sans approximation aux valeurs finales (P1 6,67 bar ; P2 3,33 bar), ici égales aux valeurs initiales. En prenant des valeurs de P1 et P2 nulles comme conditions initiales, le premier pas d'itération à l'instant t1 donne les valeurs de gains et de pressions suivantes à partir des relations II : K1 = K/ PA , K2 = K3 = K/sm, PN1 = I PA, P2 = P1/2 1+ 12 /2sm En généralisant à un pas de calcul n, correspondant à un instant tn, les valeurs de gains Kin, Ken et Kan, ainsi que celles des pressions Pin et Pen, s'expriment en fonction des valeurs de gains et de pressions Plant et Plant, calculées au pas précédent (n-1) à l'instant antérieur tn_1, par les relations : K1 = K/ VPA -Plant K2 = K3 = K/Plant -P2ant P1 = 1 1+jPA-Plant / 2V 11ant-P2ant Dans cet exemple, les valeurs restent quasiment stables après dix pas de calcul, avec un écart ultérieur entre deux pas inférieur à 5 %o. Afin d'appliquer le procédé de l'invention à un circuit complet de type de commande de direction - pour lequel un traitement matriciel est plus adapté pour résoudre le système d'équations linéarisées - l'invention est d'abord illustrée ci-après en figure 2 pour une portion de circuit simplifiée du circuit de commande de direction, comportant un générateur de pression 20. Un tel générateur de pression intervient dans les circuits hydromécaniques, en particulier dans le circuit de commande de direction qui sera décrit ultérieurement. L'équation au noeud N3 de la portion de circuit 3 s'écrit, en appliquant la linéarisation par corde dans les relations des débits Q1 et Q2 en fonction des pressions PS, PR et P3 (PS étant la pression fournie par le générateur de pression 20, PR la pression de retour à la bâche 22, et P3 la pression au noeud N3) et des gains KS et KR des restricteurs 21 et 23, respectivement associés au générateur et au retour à la bâche : PA, P2 = Plana KS (Ps-P3)- KR (P3-PR) (IV) VPS -P3ant VP3ant -PR En introduisant, comme précédemment, les gains de linéarisation des restricteurs 21 et 23 : KS* (égal à KS/VPS-P3ant ) et KR* (égal à KR/VP3ant-PR ) - avec une valeur seuil de substitution en cas d'égalité des pressions au dénominateur - l'équation IV fournit l'expression de la pression au noeud P3 suivante : PSKs+PRKÂ P3 KsS+KR R Chaque pas de calcul d'itération fournit, successivement, des valeurs de gains de Iinéarisation KS*, KR* et de pression P3 par intégration de la valeur de la pression au pas de calcul antérieur - P3ant - dans le calcul des gains de Iinéarisation KS*, KR*' puis par introduction de la valeur de ces gains dans le calcul de la pression P3 à l'étape de calcul ultérieure. En prenant par exemple comme conditions initiales : PR nulle, PS égale à 10 bar, et une valeur de seuil de 10-2, l'itération fournit - après sept pas - une valeur de P3 qui excède la valeur d'équilibre (PS/2) et - après dix pas - une valeur qui diffère de la valeur d'équilibre de 9 %o. La détermination des pressions (plus généralement des potentiels) et des débits (plus généralement des flux) - à l'équilibre dans un circuit donné - est traitée par itération dans une unité de calcul. Les solutions précédentes établies pour des portions de circuits 1 à 3 sont établies à partir d'un même système d'équations linéarisées {0 faisant intervenir des gains de Iinéarisation K* de composants élémentaires. Après normalisation des potentiels Pi et flux Q; par introduction de la notation p; = Pif/Po et q; = Q;/Qo (Po et Qo étant des valeurs de référence) - afin d'améliorer la précision du calcul et la stabilité numérique -, un tel système s'écrit selon l'expression générale : {q; = E K*Op + bi} avec K* = K Les expressions de q; fournissent alors par développement des variations Op une écriture matricielle de la forme : [gHA];.[pMB], (1) q et p ayant respectivement pour coordonnées q; et p; [Al étant la matrice de modélisation d'un composant « j », et [BI un vecteur complémentaire pour le composant « j » de coordonnées b; La matrice [Al et le vecteur [BI sont respectivement une matrice carrée et un vecteur de dimensions égales au nombre de ports du composant « j ». Ils font intervenir des caractéristiques du composant visé : gain, variations de paramètres concernant une étape antérieure, dimensions, etc. Puis, en appliquant la loi des noeuds en chaque noeud « n » (Eq = 0), la relation (1) fournit une relation de sommation E ([A]i [p] + [BI) = 0, sommant en chaque noeud les contributions de tous les composants connectés par ce noeud. Par intégration matricielle, cette sommation s'écrit alors, en chaque noeud « n » : [Ain [p] + [B]n = 0 (2) Les matrices [Ain et vecteurs [B]n de contribution des composants du circuit en chaque noeud « n » sont intégrées - par pavé - dans des matrices globales de circuit [A]c et [B]c, selon une topologie de position des noeuds propre au circuit. La relation (2) fournit alors les valeurs de potentiels par la résolution matricielle : ( [p] = - [B]/[A])c (3) La position des noeuds est définie pour chaque composant en fonction des ports par une matrice de topologie [N]c. Les valeurs de potentiel aux ports d'un composant «j » donné sont alors fournies par extraction des potentiels [pic de la relation (3) correspondant aux noeuds identifiés dans la matrice de topologie [N]c. Les valeurs de flux sont alors fournies par la relation (1). La figure 3 illustre l'exemple d'application du procédé de l'invention sur un circuit 4 de commande de direction hydromécanique de roulette d'avion, appelé « circuit de commande roulette » ou « circuit CR ». Ce circuit comporte des composants hydromécaniques couplés entre eux de la façon suivante : - un générateur de pression de fluide avec retour à la bâche, 30, du type décrit en référence au circuit élémentaire 3 de la figure 2, - une servovalve 31 de répartition du fluide entre - deux valves anti-shimmy 32a et 32b, comportant chacune une valve de contrôle 33 et un restricteur 34 montés en parallèle, chaque valve 32a et 32b étant couplée respectivement à l'une des deux chambres de compression à capacité variable 35a et 35b - un vérin 35 d'un actionneur 36, qui comporte un piston 35c définissant la capacité des chambres 35a et 35b et couplé à un pignon 37 sur crémaillère 37c, le pignon 37 étant en liaison avec - une colonne de direction 38 de roulette, composée d'une roue d'inertie 38a couplée à un générateur de couple antagoniste 38b de réaction au sol : l'inertie augmente par déplacement axial de la roue 38a et le générateur antagoniste se compose, dans l'exemple illustré, d'un amortisseur rotatif 39a intégrant un ressort hélicoïdal 39b. Un asservissement de la position de la servovalve 31 est réalisé par une boucle rétroactive 3b de la mesure de l'angle 0 de direction. La mesure est effectuée par un capteur 40 sur la colonne 38 en comparaison 41 avec une valeur de consigne 42. Selon le procédé de l'invention, une modélisation de chaque composant d'un circuit hydromécanique donné est appelée d'une librairie de composants apte à fournir les matrices de modélisation [Al et [BI de chaque composant «j » du circuit. Ainsi, pour chaque composant du circuit 4, les matrices fournies par la librairie définissent les expressions combinées des caractéristiques du composant à saisir : gains intrinsèques et gains de linéarisation pour tous les composants, ouverture pour la servovalve 31, pression de craquage pour les valves anti-shimmy 32a et 32b, dimensions du vérin 35 et du pignon 37 pour l'actionneur 36, ainsi qu'inertie en rotation et coefficient d'amortissement pour la colonne de direction 38. Comme paramètres de flux et de potentiel, le débit et la pression interviennent pour tous les composants sauf la colonne de direction. Pour l'actionneur 36, par son pignon 37, et la colonne de direction 38 font intervenir le couple exercé sur l'élément en rotation ainsi que sa vitesse angulaire, respectivement comme paramètres de type flux et potentiel. Ainsi : - pour le générateur de pression 30, selon les notations déjà utilisées pour ce composant : Ksx 0 LAIg = 0 KR ' - pour la servovalve 31: LBIg - Ksxps - KR,,pR ~K13+K14 o -1(13 - K14 * K* 24 x o Kz3+K* -K23 K 24 5 SV - - K13 Kz3 K13+K _ 23 , x K14 K24 o K14+ K2 (les indices renvoient à la numérotation des ports de la servovalve) 20 - pour les valves anti-shimmy 32a et 32b dites « as » en abrégé, avec une valve de contrôle de gain Kcv, si 1) la variation de pression Op entre les ports du composant est supérieure ou égale à la pression de craquage pc de la valve de contrôle : 25 [A]as = (Kcv + KR) -1 -1 Kcv pc 1 [B]as= - Kcv pc 2) dans le cas contraire (Op 30 [A]as= KR[ 1 1 ' [B]as = 0 - pour l'actionneur 36, muni du vérin 35 de sections internes s1 et s2 des chambres 35a et 35b, et du pignon 37 de rayon « r », de couple « c » et de vitesse de rotation angulaire « co », la relation (1) [gHA]i.[p]+[B]i a comme éléments : 0 0 rs1 0 0 - rs2 ~rs1 - rs2 0 J q1 q2 C , [AL= [G = [p]a= , [B]a = 0 - pour la colonne de direction 38, comportant un générateur de couple antagoniste 38b de coefficient d'amortissement « 1 » et de constante de raideur « k », la relation (1) s'écrit sous une forme linéaire c = Acrco + Bcr avec : Acr = - J/dt -1 - kdt, Bcr = Jutant/dt - niant, gant étant la mesure de l'angle de rotation à l'étape de calcul antérieur et J l'inertie en rotation. Un exemple de construction de matrice de topologie de circuit est maintenant décrit par une numérotation élémentaire des ports, des composants et des noeuds de ces composants puis positionnement dans la matrice. Un circuit hydromécanique modélisé comporte « n » noeuds identifiés par une numérotation de 1 à n. Les ports d'entrée/sortie sont identifiés par composant selon une numérotation de 1 à k, et les composants sont classés selon une numérotation 1 à j. La matrice de topologie [N] est alors construite pour représenter l'architecture du circuit hydromécanique. Cette matrice a un nombre de lignes égal au nombre de composants et un nombre colonnes égal au nombre de ports du composant ayant le plus grand nombre de ports. Les éléments de cette matrice sont constitués par les numéros des noeuds, chaque numéro de noeud étant positionné avec pour coordonnées : en ligne, le numéro du composant du noeud et, en colonne, le numéro du port auquel le noeud est connecté. Ainsi, dans le cas du circuit 4, sept noeuds n1 à n7 sont identifiés entre les différents composants et sont numérotés élémentairement 1 à 7. Les ports de chaque composant sont identifiés par une numérotation élémentaire indépendante, à savoir : 1 et 2 pour les ports L1 et L2 du générateur de pression 30, 1 à 4 pour les ports L3 à L6 de la servovalve 31, 1 et 2 pour les ports L7, L8 et L9, L10, respectivement des valves anti-shimmy 32a et 32b, 1 à 3 pour les ports L11 à L13 de l'actionneur 36, et 1 pour le port d'entrée L14 de la colonne de direction 38. Quant aux six composants - respectivement : le générateur 30, le servovalve 31, les valves anti-shimmy 32a et 32b, l'actionneur 36 et la colonne de direction 38 - ils sont numérotés élémentairement 1 à 6. La matrice de topologie du circuit CR, [N]cR, a alors pour éléments les numéros des noeuds 1 à 7 placés selon les numéros de ports 1 à 4 en colonne et les numéros de composants 1 à 6 en ligne, tels que définis ci-dessus, ce qui fournit une matrice de la forme: 1 6 1 6 2 5 2 3 5 4 3 4 7 7 J [N]cR = Par itération à partir de la consigne d'orientation roulette et des matrices [Al et [BI fournies pour tous les composants, le calcul permet de converger quasi-instantanément vers des valeurs d'équilibre pour les éléments de ces matrices. Les matrices globales du circuit, [A]cR et [B]cR sont alors remplies au fur et à mesure de la détermination des matrices de chaque composant selon la topologie définie par la matrice [N]cR. Enfin, les valeurs de pressions aux ports de chaque composant sont fournies par extraction des potentiels [p]c, de la relation ([p]=- [B]/[AbcR correspondant aux numérotations de noeuds identifiés dans la matrice de topologie [N]cR pour le composant et pour chaque port de ce composant. Par exemple, pour l'actionneur 36 (5è' composant dans la matrice de topologie), les numéros de noeuds sont 3, 4 et 7 aux ports 1, 2 et 3. Les valeurs de débits sont ensuite fournies par la relation (1). Le facteur de compressibilité du fluide dans les conduits a été négligé dans les exemples décrits ci-dessus qui expriment la loi des noeuds sous la forme ZQ = 0, car ce facteur de compressibilité reste de valeur faible par rapport aux autres caractéristiques, par exemple de gain. Afin d'améliorer en représentativité les performances du procédé, sans nuire à la compatibilité au temps réel, il peut être avantageux de faire intervenir ce facteur de compressibilité du fluide dans la loi des noeuds, et de le linéariser selon une approche par schéma implicite telle que présentée ci-dessous. L'expression de la sommation des débits gQ) aux noeuds en fonction des variations de pression a été linéarisée dans les exemples précédents sous la forme : E Q(P) = K*OP + B L'introduction du facteur de compressibilité du fluide, qui fait intervenir la dérivée de sa pression « P » en fonction du temps « t », permet d'exprimer la loi des noeuds sous sa forme rigoureuse, à savoir (avec la convention de signe positif pour les débits sortants) : ZQ(P) = - R dt ' avec V = volume du conduit, p = module de masse La linéarisation de l'équation aux noeuds, après introduction du facteur de compressibilité, est mise en oeuvre par une approche tangentielle qui substitue la dérivée partielle de potentiel en fonction du temps (dP/dt) par une variation de ce potentiel en fonction de la variation de temps finie correspondante. Cette approche est conduite ici selon le schéma implicite, qui permet d'exprimer l'approximation de la dérivée par : dP P-P* dt At Après introduction de cette approximation et développement, l'équation aux noeuds avec schéma implicite se présente sous la forme d'une fonction linéaire de la pression f(P) constante (ci-après « Cste ») : zQ(P) + p.& P R Ot P Ou f(P) = Cste Le choix du schéma implicite pour formaliser l'approximation tangentielle permet de converger dans la détermination des pressions et débits en faisant un choix de traitement plus complexe - car il introduit des facteurs supplémentaires (P*) - mais totalement gérable par l'unité numérique, tout en apportant plus de précision et assurant la stabilité de la résolution. La prise en compte de la compressibilité du fluide permet d'intégrer un facteur d'amortissement réaliste et ainsi de rendre compte d'adaptations transitoires qui sinon pourraient se traduire en simulation, dans de rares cas, par des situations de blocage. Par exemple dans la commande de colonne de direction de roulette de nez de gros aéronefs (systèmes dits de « push-pull » Ou « pousser-tirer »), les axes des bras d'orientation qui sont des vérins peuvent passer par l'axe de rotation de la colonne. Des situations de blocage peuvent apparaître en simulation car aucun couple ne peut alors s'exercer sur la colonne. La prise en compte de la compressibilité du fluide permet d'éliminer ces états de blocage en simulation. Ainsi, en ce qui concerne l'actionneur 36, la prise en compte du coefficient de compression permet d'écrire la matrice de modélisation [A], et le vecteur complémentaire [I3], sous la forme : KR 0 rs i - KR [A], = 0 KR -rs2 , [B]a = - KR Rsi - rs2 0 0 La présente invention n'est pas limitée aux exemples décrits et représentés. Sa portée s'étend à tout système physique de complexité pouvant atteindre un haut degré, tout système physique étant décomposable en portions élémentaires d'un point de vue de son comportement simulé conformément à la modélisation selon l'invention. REVENDICATIONS1. Procédé de simulation de comportement d'un système physique complexe en temps réel dans une unité de calcul numérique par des étapes de simulation stables de durée suffisamment représentative des phénomènes physiques dans le domaine fréquentiel, caractérisé en ce qu'il consiste à créer une librairie de composants, à assembler des composants de cette librairie comportant au moins un composant de comportement non linéaire pour modéliser une architecture de circuit par interconnexion des ports d'entrée des composants à travers des noeuds d'accès, chaque port étant connecté à un seul noeud, à caractériser chaque port de chaque composant par un couple de - paramètres potentiel/flux, à linéariser le comportement du (des) composant(s) non linéaire(s) selon une approximation par corde en un nombre de pas réduit permettant une convergence dans un laps de temps réduit, et à établir une relation de détermination du flux à partir de la-variation de potentiel à chaque port des composants de comportements linéaires et linéarisé(s), le potentiel de chaque port étant assimilé à celui de son noeud de connexion, à constituer un système d'équations linéaires et linéarisée(s) entre les paramètres en appliquant la loi des noeuds en chaque noeud et à résoudre ce système d'équations de manière itérative pour déterminer les valeurs des paramètres à chaque noeud. 2. Procédé de simulation selon la revendication 1, dans lequel le circuit comporte au moins un composant non linéaire intemporel, dont le comportement n'est pas lié par une variation des paramètres en fonction du temps, en particulier par une dérivée de paramètre en fonction du temps. 3. Procédé de simulation selon l'une des revendications 1 ou 2, dans lequel les paramètres de potentiel sont choisis parmi : la pression, la vitesse angulaire et la vitesse linéaire, et les paramètres de flux choisis parmi le débit, le couple et la force. 4. Procédé de simulation selon l'une quelconque des revendications précédentes, dans lequel les composants non linéairesintemporels sont choisis parmi des générateurs de pression, des servovalves, des valves, des restricteurs de flux, des vérins et des colonnes de direction de véhicule aérien ou terrestre. 5. Procédé de simulation selon l'une quelconque des revendications précédentes, dans lequel la loi des noeuds fait intervenir un facteur de compressibilité temporel à linéariser selon une approche par linéarisation tangentielle avec schéma implicite. 6. Procédé de simulation selon l'une quelconque des revendications précédentes, dans lequel le système d'équations linéaires et linéarisée(s) entre le paramètre de flux et le paramètre de potentiel pour chaque noeud autour d'un composant Cj, à partir des variations du paramètre de potentiel à un instant antérieur selon l'approximation par corde pour le(s) comportement(s) non linéaire(s), est transposé en mode matriciel par formation d'une matrice de modélisation [A]; et d'un vecteur complémentaire pl, puis pour l'ensemble des composants Cj du circuit par intégration des matrices [A]j dans une matrice de modélisation globale [A] et des vecteurs complémentaires [B], dans un vecteur global [B] selon une organisation définie par une matrice de topologie [N] du circuit. 7. Procédé de simulation selon l'une quelconque des revendications précédentes, dans lequel la topologie d'un circuit est modélisée en numérotant chaque composant, chaque port de chaque composant et chaque noeud, puis en positionnant les numérotations des noeuds en mode matriciel dans la matrice [N] de topologie en ligne selon la numérotation des composants et en colonne selon la numérotation des ports. 8. Procédé de simulation selon l'une quelconque des revendications précédentes, dans lequel - un circuit étant modélisé - la simulation de son comportement par les valeurs de flux aux ports de chaque composant Cj à partir des valeurs de potentiel à chaque noeud par application de la loi des noeuds, ést déterminée par des matrices [Al, et vecteurs [B], de contribution des composants du circuit en chaque noeud « n » sont intégrées - par pavé - dans des matrices globales de circuit [A]c et [Mc, selon la topologie de position des noeuds propre au circuit, les valeurs de potentiels étant fournies par la résolution matricielle ([p]=-[BMA])c, et la position des noeuds étant définie pour chaque composant en fonction des ports par la matrice de topologie [Mc, les valeurs de potentiel aux ports d'un composant « j » donné sont alors fournies par extraction des potentiels [pi, correspondant aux noeuds identifiés dans la matrice de topologie [N]c, les valeurs de flux sont alors fournies par la relation flux/potentiel.