Brevet FR3145412A1 — publié le 20240802

Procédé de prédiction de la propagation d’une fissure par fatigue dans une pièce mécanique comprenant : une étape de maillage à l’instant t 0 ; une étape de détermination (802) d’un champ de contraintes internes (9032)) ; une réduction de modèle (803) ; un calcul de facteurs d’intensité des contraintes modaux (806) le long d’un front de fissure ; une reconstitution (807) de l’évolution temporelle des FIC ; une estimation (808) d’une avancée de fissure ; l’obtention d’un maillage fissuré propagé ; les étapes b. à g. étant répétées pour chaque pas de la propagation modélisée, de l’instant t 0 + Δ t jusqu’à un instant final t f . Figure pour l’abrégé : Fig. 5 Procédé de calcul de facteurs d’intensité des contraintes et de prévision d’une propagation de fissure dans une pièce mécanique DOMAINE DE L'INVENTION L’invention concerne de manière générale le domaine de l’analyse de la propagation de fissures dans des pièces mécaniques. L’invention se situe dans le domaine des méthodes de simulation. Elle vise à estimer la tolérance aux dommages des pièces mécaniques, ainsi qu’à calculer des durées de vies associées à la propagation de fissures, dans le cadre d’une propagation par fatigue avec une série de cycles de chargement. Plus particulièrement, l’invention concerne des méthodes de détermination et d’interpolation de facteurs d’intensité des contraintes (abrévié par la suite en FIC) de pièces mécaniques modélisées numériquement. Elle trouve très avantageusement application dans le domaine aéronautique, notamment pour le suivi de pièces de turbomachines fortement sollicitées mécaniquement et thermiquement, mais peut également être utilisée dans d’autres domaines. Dans le domaine aéronautique, on cherche notamment à pouvoir estimer la durée de vie restante d’une pièce à compter de la détection d’une anomalie dans cette pièce, notamment en termes de nombre de missions en vol que la pièce pourra effectuer dans le futur. L’invention concerne en outre l’application de telles méthodes dans le cadre de la maintenance et de la réparation aéronautique, afin de permettre la prévision de la mise hors service ou du remplacement de pièces en fonction des durées de vie obtenues. Elle concerne également la démonstration, par calcul, de la tenue en propagation de la pièce modélisée, en vue de sa certification et/ou de l’extension de sa durée de vie. ETAT DE LA TECHNIQUE Dans le domaine de la mécanique de la rupture, il est connu de décrire l’intensité des champs mécaniques asymptotique au voisinage de l’ensemble des points décrivant le front de fissure à l’extrémité d’une fissure grâce à des facteurs d’intensité des contraintes (FIC dans la suite du texte) au sein de la pièce qu’on souhaite analyser. Pour une dimension de fissure donnée et en un point donné du front de fissure, le FIC est la grandeur caractérisant l’intensité du chargement au voisinage de la pointe de fissure. Les FIC se composent de trois grandeurs, notées KI, KII et KIII, correspondant respectivement aux trois modes d’ouverture de la fissure, de cisaillement plan et de cisaillement hors plan. Selon l’amplitude de ces grandeurs au cours d’un cycle de chargement de fatigue, la fissure va se propager plus ou moins vite et la forme du front de fissure va évoluer. Le rôle des simulations numériques de propagation est de fournir les valeurs des FIC le long du front de fissure pour différents pas de propagation, et donc pour différentes géométries de fissure. Dans l’état de la technique, on distingue deux catégories d’approches de calcul des durées de vie de pièces mécaniques en propagation de fissure : les approches dites « fissuration 2D » et les approches dites « fissuration 3D ». Méthodes fissuration 2D En fissuration 2D, on déduit le champ de contraintes de la structure fissurée par post-traitement semi-analytique du champ de contrainte de la structure non-fissurée, qui doit donc être une donnée du problème. Pour ces méthodes, le champ de contraintes σ(x,y) est obtenu en sélectionnant un nœud du maillage pour représenter le centre de la fissure, en orientant le plan de propagation correctement et en adaptant ses dimensions de façon à être représentatif de la structure. Un tel post-traitement consiste à interpoler le champ de contraintes au niveau d’un plan de propagation qui traverse la structure non-fissurée. Ce plan de propagation, restreint sur une zone d’intérêt (appelée barreau ou plan de propagation), a des dimensions telles qu’elles sont représentatives de la géométrie réelle de la structure. Le plan de propagation est orienté de manière à ce que la direction de la contrainte principale maximale, c’est à dire maximisée sur l’intervalle de temps considéré, soit orthogonale à ce dernier. Par « contrainte principale maximale », on entend la valeur propre la plus élevée du tenseur de contraintes, cette contrainte étant dirigée selon une direction orthogonale au plan tangent à la fissure et passant par son centre. Le champ de contrainte sur le plan de propagation est ensuite post-traité pour obtenir le FIC correspondant à un instant donné de la mission et pour une géométrie donnée de fissure. Il existe principalement deux types de post-traitements permettant le calcul des FIC. Premièrement, la méthode dite « polynomiale » consiste en une discrétisation du plan de propagation en une grille 2D, par exemple de 100x100 points. Pour chaque instant de la mission, on identifie un polynôme de contrainte grâce à l’ensemble des contraintes associées à chacun des nœuds de la grille. Ainsi, pour toute évolution temporelle des contraintes sur la grille, le polynôme de contrainte ainsi défini sera modifié. Pour un champ de contrainte unidimensionnel, le polynôme de contrainte s’écrit où x est la position, variant entre 0 et la profondeur du plan de propagation T, t est l’instant de la mission considéré et a i sont les coefficients du polynôme. Le facteur d’intensité des contraintes peut alors se déduire par linéarité du problème où geom correspond à l’ensemble des paramètres donnant la géométrie de la fissure ainsi que du plan de propagation, et K I i (geom,t) est la valeur du FIC pour un champ de contrainte appliqué sur le plan de propagation a i (t).(x/T) i . Ces paramètres sont notamment détaillés dans la référence [6]. Une deuxième méthode de calcul des FIC, dite méthode des « fonctions de poids », utilise une « fonction de poids » m qui dépend de la position dans le champ de contrainte ainsi que des paramètres de géométrie de la fissure et du plan de propagation. On calcule l’intégrale sur la surface concernée : Pour le secteur aéronautique, les méthodes de fissuration 2D présentent l’avantage d’être faciles à mettre en œuvre et de permettre des calculs rapides, une simulation de mission en vol pouvant durer de quelques secondes à quelques minutes. Cependant, elles présentent plusieurs inconvénients : elles ne permettent d’effectuer des simulations que pour des géométries simples, et ne permettent pas de traiter des coalescences de fissures (c’est-à-dire des cas où plusieurs fissures interagissent entre elles). Dans le cas où les contraintes ont une orientation qui change au cours du temps, l’hypothèse de propagation plane qui est faite dans ces méthodes n’est plus réaliste. Les méthodes de fissuration 2D ne tiennent compte que du mode I (mode d’ouverture de la fissure), et ne permettent pas d’obtenir des résultats fiables lorsque le chargement réel est de mode mixte, c’est-à-dire comprenant des composantes selon les modes II et III. Enfin, un inconvénient majeur des méthodes de fissuration 2D est qu’elles sont trop conservatives, et peuvent mener à des résultats de durée de vie qui ne satisfont pas aux objectifs prévus, non pas car la structure réelle romprait dans des conditions de chargement similaires mais à cause des hypothèses conservatives utilisées pour la modélisation. Ce conservatisme est notamment dû au fait qu’on impose toujours une contrainte en tant que conditions limites, ce qui ne permet notamment pas de prendre en compte l’assouplissement de la structure au cours de la propagation de la fissure, par exemple lorsque la structure est soumise à des déplacements imposés plutôt qu’à des contraintes imposées, cet assouplissement réduisant les efforts vus par la structure. La illustre le fonctionnement des procédés de fissuration 2D. Les données d’entrée 100 d’une telle simulation sont les données relatives à la fissure 101, le champ de contraintes évolutif sur le plan de propagation 102, le champ de température 103 et la loi de propagation de la fissure 104. L’algorithme 200 de calcul de la durée de vie permet d’obtenir les résultats 300, à savoir le nombre de cycles effectués jusqu’à la rupture 301 et les données relatives à la fissure finale 302. Les résultats de post-traitement donnent l’évolution morphologique de la fissure au cours de la propagation. Une boucle débute par la détermination des FIC à l’étape 202 pour chaque instant du chargement et pour chaque point du front de fissure, au moyen de l’équation 3. Puis, à l’aide d’un traitement du chargement 203 associé à la loi de propagation, il est possible d’en déduire l’incrément d’avancée de fissure 204 pour chaque point du front de fissure. Ensuite, on vérifie qu’aucun critère de rupture 205 n’est atteint. Si un critère de rupture est atteint, l’algorithme s’arrête et le nombre de cycles 301 réalisés jusque-là correspond à la durée de vie en propagation de fissure. Sinon, l’algorithme 200 se poursuit. On effectue alors une mise à jour 206 des dimensions de fissure avec les incréments d’avancées de fissure 204 calculés précédemment. La nouvelle géométrie obtenue servira pour les calculs du cycle suivant. Avant de débuter le cycle suivant, il convient néanmoins de vérifier les critères de transition, à l’étape 201. En effet, il est possible que la fissure quitte le périmètre d’applicabilité d’une équation du type représentée par l’équation 3, c’est par exemple le cas d’une fissure proche d’un bord. Pour que la propagation puisse se poursuivre, il suffit alors de changer la géométrie de fissure ainsi que l’équation associée. Une telle transition se fait de manière conservative, et s’accompagne d’une augmentation de la surface fissurée. Ceci conduit à des prédictions de durée de vie de pièces conservatives par rapport à la réalité. Méthodes fissuration 3D En référence à la , dans les méthodes 3D, les données d’entrée 400 comprennent une géométrie de la fissure initiale 401, le champ de contraintes sur la structure saine ou non-fissurée 402, la loi de propagation utilisée 403 et le champ de température 404 auquel est soumise la pièce. Ensemble, les champs de contrainte 402 et de température 404 constituent les conditions aux limites du problème à résoudre. Le procédé de propagation 500 de la fissure commence par une étape de maillage 501 de la structure fissurée. Le maillage peut être fait de manière explicite, par exemple par la méthode des éléments finis, ou de manière implicite, par exemple par la méthode des éléments finis étendus. Après l’étape de maillage 501, on utilise une méthode de résolution numérique appropriée dans une étape d’obtention du FIC 502. On incrémente ensuite à l’étape 503 la fissure, et on vérifie à l’étape 504 si un critère d’arrêt est atteint, auquel cas la modélisation doit s’arrêter. Si ce n’est pas le cas, on met à jour les dimensions de la fissure à l’étape 505 et on réitère le procédé à partir de l’étape 501 de maillage de la structure fissurée. Pour la résolution numérique du problème, on pourra notamment utiliser les méthodes suivantes : La méthode des éléments finis (en anglais Finite Element Method ou FEM), qui consiste à mailler explicitement la structure, y compris la fissure ; La méthode des éléments finis étendus ( Extended Finite Element Method ou X-FEM), qui consiste à mailler explicitement la structure, en maillant implicitement la fissure par le biais de fonctions de forme enrichies représentant la singularité des contraintes en pointe de fissure et la discontinuité du champ de déplacement au niveau des lèvres de la fissure. La méthode des éléments de frontière ( Boundary Element Method ou BEM), pour laquelle, pour une structure tridimensionnelle, on maille la surface de la structure plutôt que son volume. Les lèvres de la fissure sont également maillées de manière surfacique. Pour chaque géométrie de fissure, les procédés de fissuration 3D comprennent autant de résolutions numériques à effectuer que d’instants modélisés de la mission. Les méthodes 3D comprennent ensuite des étapes 600 d’obtention de la durée de vie de la pièce mécanique : Une étape de vérification 601 de l’atteinte de critères d’arrêt - atteinte d’un FIC critique KI c appelé « ténacité » ou d’une propagation de fissure maximale donnée a max ; une étape 602 d’établissement de tableaux de correspondance entre les FIC et les dimensions de fissure, puis une étape 603 de calcul de la durée de vie. Les résultats 700 fournis par une telle méthode 3D sont la géométrie finale 701 de la fissure, des informations intermédiaires 702 tels que les critères d’arrêt, la vitesse de propagation de la fissure et l’analyse des cycles, et le nombre de cycles 703. A partir des résultats concernant les facteurs d’intensité des contraintes, obtenus pour un instant de la mission (instant associé à la contrainte principale maximisée sur l’ensemble des instants, dit « instant critique »), on extrait, à l’étape 602, un tableau donnant un facteur d’intensité des contraintes équivalent en fonction d’une dimension de fissure équivalente. Il existe différentes manières de construire ce tableau. L’une d’entre elles, consiste à considérer, pour chaque front de fissure, la plus grande valeur du facteur d’intensité des contraintes. Pour la dimension de fissure équivalente, on peut considérer la plus petite distance au front du nœud représentant le centre de la fissure. L’évolution temporelle du facteur d’intensité des contraintes peut s’obtenir par proportionnalité à la contrainte principale maximale au niveau de l’anomalie. Les tableaux donnant le facteur d’intensité des contraintes en fonction de la dimension de la fissure, appelés « formulaires », peuvent alors être utilisés à la place de l’équation 3 intervenant en « fissuration 2D » pour obtenir une durée de vie. Les méthodes de fissuration 3D modélisent l’intégralité de la structure, et rendent donc compte des effets de structure. De plus, elles rendent compte de conditions aux limites plus variées : on n’impose pas nécessairement une contrainte, mais on peut prévoir des déplacements ou forces imposés, des encastrements, contacts ou champs thermiques. Les méthodes 3D présentent également l’avantage de pouvoir modéliser des fissures non planes, de géométries parfois complexes, car ces méthodes numériques font intervenir explicitement la fissure dans la modélisation. Enfin, les méthodes 3D prennent en compte les trois modes de propagation de la fissure, alors que les méthodes 2D se basent sur l’hypothèse que le chargement est intégralement en mode I. L’inconvénient majeur des méthodes 3D est la complexité de la mise en œuvre, avec des temps de calculs beaucoup plus longs, pouvant aller jusqu’à plusieurs mois pour la simulation complète d’une mission en vol. De plus, les fichiers de résultats sont plus lourds et leur post-traitement est plus complexe. En outre, les méthodes 3D fonctionnent sur le principe d’une mise à jour du maillage de la pièce à chaque étape de propagation, et il est connu que ceci rend difficile l’utilisation d’une réduction de modèle de type modal - par exemple, par décomposition orthogonale propre - qui permettrait de réduire le nombre de calculs à effectuer pour obtenir les FIC. Autres méthodes Outre les méthodes de fissuration 2D et 3D, il existe d’autres méthodes permettant le calcul du facteur d’intensité des contraintes. On peut citer la méthode dite des « dislocations distribuées », qui s’apparente à une méthode de type éléments de frontières, mais les fonctions de formes utilisées pour résoudre le problème mécanique sont inspirées des champs de contraintes induits par la présence de dislocations dans un milieu linéaire, homogène et isotrope. Cette méthode permet le calcul très rapide des facteurs d’intensité des contraintes en modes mixtes et pour des géométries de fissures très complexes. Cependant, elle est limitée aux milieux infinis et au milieux semi-infinis. Un but de l’invention est de fournir un procédé automatique de calcul des FIC, sans avoir à intervenir dans le code de calcul afin de garantir une mise en œuvre rapide. Un autre but de l’invention est d’obtenir une chaîne complète de calcul numérique simplifiée d’une durée de vie en propagation de fissure pour une pièce métallique sous un chargement de fatigue, notamment dans le cas d’une propagation par cycles de missions. Un autre but de l’invention est de permettre d’étendre la durée de vie mécanique des pièces, en permettant de savoir avec précision quand une pièce doit être changée. Un autre but de l’invention est de fournir un procédé d’estimation de la durée de vie d’une pièce plus fiable et moins conservatif que les procédés de fissuration 2D et permettant la prise en compte de missions complexes, éventuellement avec des chargements non-proportionnels (c’est-à-dire dont les gradients de contrainte normalisés évoluent avec le temps et dont les directions principales évoluent avec le temps et l’espace) tout en étant plus facile à mettre en œuvre et moins cher que les procédés de fissuration 3D. Ceci permettra notamment la réalisation de simulations en bureau d’étude. Il est à cet effet proposé, selon un premier aspect de l’invention un procédé de prédiction de la propagation d’une fissure par fatigue dans une pièce mécanique soumise à un champ de contraintes externe et modélisée numériquement entre un instant initial t 0 et un instant final t f avec un pas de Δt, comprenant les étapes suivantes : une étape de maillage à l’instant t 0 ; une étape de détermination d’un champ de contraintes internes à appliquer aux lèvres de la fissure, par calcul du champ de contraintes aux mêmes points du maillage pour la structure non fissurée lorsque le champ de contraintes externes lui est appliqué, le champ de contrainte interne de la pièce réelle fissurée étant alors égal et de signe opposé au champ de contrainte interne de la pièce non fissurée équivalente ; caractérisé en ce que l’étape de maillage comprend les sous-étapes suivantes : Maillage d’une pièce non fissurée identique à la pièce fissurée ; Maillage numérique de la fissure ; Insertion du maillage de fissure dans le maillage de la pièce non fissurée, de sorte à obtenir un maillage numérique de la pièce fissurée ; le procédé comprenant en outre une étape de réduction de modèle de la structure réelle à un modèle modal simplifié, la réduction de modèle consistant à exprimer le champ de contrainte à chaque nœud du maillage fissuré sous forme d’une combinaison linéaire de champs de contraintes modaux (σ k ) associés chacun respectivement à un mode structurel k de la pièce, pour un nombre de modes prédéterminé ; une étape de calcul de facteurs d’intensité des contraintes modaux le long d’un front de fissure pour le modèle modal simplifié ; une étape de reconstitution de l’évolution temporelle des FIC à chaque nœud du front de la fissure à partir des FIC modaux obtenus à l’étape de calcul ; une étape d’estimation (808) de la valeur d’une avancée de fissure à partir des évolutions temporelles des FIC obtenues à l’étape de reconstitution et au moyen d’une loi de propagation de fissure, l’étape d’estimation comprenant la détermination de la direction d’avancée de chaque nœud du front de fissure par calcul de l’angle qui maximise la vitesse de propagation de ce nœud ; une étape d’obtention d’un maillage fissuré propagé reconstitué à partir du maillage de la pièce fissurée et de l’estimation de la valeur d’avancée de fissure ; les étapes b. à g. étant répétées pour chaque pas de la propagation modélisée, de l’instant t 0 + Δ t jusqu’à l’instant final t f . Un tel procédé est en outre avantageusement complété par les différentes caractéristiques suivantes prises seules ou selon toutes leurs combinaisons techniquement possibles : - le procédé comprend en outre : une étape d’établissement de formulaires consistant à établir, pour chaque pas de propagation, des tableaux de correspondance entre des facteurs d’intensité de contraintes calculés et des dimensions de fissure ; une étape de calcul d’une durée de vie restante de la pièce mécanique sur la base des formulaires établis ; - la réduction de modèle est une décomposition orthogonale propre (POD) ; - la loi de propagation de fissure utilisée dans l’estimation de la valeur d’avancée de fissure étant une loi de Paris ; - la loi de propagation de fissure utilisée dans l’estimation de la valeur d’avancée de fissure étant une loi de Paris à laquelle une correction est appliquée ; - l’étape d’estimation de la valeur d’avancée de fissure comprenant l’imposition, à chaque pas de propagation, d’une valeur de propagation maximale a max de la fissure de façon à réduire le nombre de pas de propagation calculés pour cette étape. - le maillage de la structure fissurée étant décomposé en : une zone remaillée, à proximité de la fissure, pour laquelle les étapes d’obtention d’un maillage numérique de la fissure et d’obtention d’un maillage numérique de la pièce fissurée sont renouvelées à chaque pas ; une zone non-remaillée, en dehors de la zone remaillée, pour laquelle le maillage est déterminé à l’étape d’obtention d’un maillage numérique de la pièce non fissurée et est maintenu pour l’ensemble des pas de propagation de la fissure ; et les étapes de calcul de facteurs d’intensité de contrainte le long du front de fissure pour chaque mode structurel du modèle simplifié, d’obtention d’un maillage fissuré propagé et le cas échéant de reconstitution du facteur d’intensité des contraintes à tous les instants de la modélisation ne sont effectuées que pour les degrés de liberté de la pièce correspondant à la zone remaillée ; L’invention concerne en outre un procédé de maintenance d’une pièce d’aéronef entre une première série (Mn) de missions en vol et une deuxième série (Mn+1) de missions en vol, comprenant les étapes suivantes : - Inspection de la pièce afin de déterminer si une anomalie susceptible de causer l’apparition d’une fissure est apparue dans la pièce suite à la première série (M n ) de missions en vol Si oui, mise en œuvre du procédé précité afin de déterminer une durée de vie de la pièce ; Sinon, remise en service de l’aéronef pour la deuxième série de missions en vol ; L’étape de détermination d’une durée de vie peut être suivie d’une comparaison avec une durée de vie de la pièce non fissurée et/ou de comparaison avec une durée de vie minimum nécessaire à une certification de la pièce. L’invention concerne également un procédé d’établissement d’une base de données pour une pièce mécanique, reliant des données d’entrée à des données de sortie ; les données d’entrée comprenant des tailles et/ou géométries de fissures, des positionnements de fissures dans la pièce et des contraintes externes à la pièce ; les données de sortie comprenant des facteurs d’intensité des contraintes et/ou des durées de vie de la pièce, comprenant les étapes suivantes : Mise en œuvre du procédé de prédiction de la propagation d’une fissure par fatigue dans une pièce mécanique tel que défini précédemment, de sorte à obtenir des facteurs d’intensité des contraintes et/ou une durée de vie de la pièce mécanique pour les conditions d’entrée ; Enregistrement, dans la base de données, d’une correspondance entre les données d’entrée et les données de sortie. Selon un autre objet, l’invention concerne un système comprenant une unité de traitement comprenant des moyens de calcul et une mémoire, l’unité étant configurée pour mettre en œuvre un des procédés tels que définis précédemment. Elle concerne également un produit programme d’ordinateur configuré pour être mis en œuvre par un système selon la revendication précédente, et comprenant des instructions pour provoquer la mise en œuvre d’un des procédés tels que définis précédemment. DESCRIPTION DES FIGURES D’autres caractéristiques, buts et avantages de l’invention ressortiront de la description qui suit, qui est purement illustrative et non limitative, et qui doit être lue en regard des dessins annexés sur lesquels : La est un schéma fonctionnel montrant les étapes du procédé de calcul d’une durée de vie de pièce mécanique fissurée selon les méthodes « fissuration 2D » de l’état de la technique. La est un schéma fonctionnel montrant les étapes du procédé de calcul d’une durée de vie de pièce mécanique fissurée selon les méthodes « fissuration 3D » de l’état de la technique. La illustre un front de fissure pour différents pas de propagation. La illustre un système permettant de mettre en œuvre l’invention, ainsi qu’une pièce à analyser. La est un schéma fonctionnel résumant les étapes du procédé de prévision de la propagation d’une fissure selon l’invention. La illustre le principe de superposition utilisé dans le procédé selon l’invention. La illustre l’utilisation de zooms structuraux, mis en œuvre dans le cadre des méthodes « fissuration 3D » de l’état de la technique. La illustre une avancée de fissure, pour une fissure subissant une extension de fissure avec bifurcation. Les et illustrent une réduction dimensionnelle intervenant dans le procédé selon l’invention. La illustre une décomposition de la structure en une zone remaillée et une zone non-remaillée. La est un schéma fonctionnel résumant les étapes du procédé de maintenance d’une pièce selon l’invention. Sur l’ensemble des figures, les éléments similaires portent des références identiques. DESCRIPTION DETAILLEE DE L'INVENTION La illustre un front de fissure pour différents pas de propagation, avec les points représentés correspondant à différents endroits du front de fissure. En référence à la , la pièce à analyser 20 est généralement une pièce destinée à l’aéronautique, pour laquelle la durée de vie doit pouvoir être estimée. La pièce 20 est typiquement une aube de fan de compresseur ou de turbine, un disque moteur, une bride, une attache moteur, ou un carter. Il peut de manière générale s’agir de toute pièce métallique, et cette liste est donnée à titre illustratif, puisque, pour la mise en œuvre du procédé, le type de pièce n’intervient cependant pas. La représente également un système 10 d'interpolation de données concernant une pièce à analyser 20, qui est modélisée numériquement. Comme indiqué en introduction, le FIC K se décompose en trois données K l , K II , K III . Le système 10 comprend une unité de traitement de données 12, par exemple un ordinateur ou un serveur, disposant de moyens de calcul 14 configurés pour mettre en œuvre des procédés qui seront décrits plus en détail ci-après et, disposant avantageusement, d'une mémoire 16. Les moyens de calcul 14 peuvent par exemple être un calculateur de type processeur, microprocesseur, microcontrôleur, etc. La mémoire 16 peut être par exemple un disque dur, une mémoire « flash », ou un espace de stockage délocalisé, de type « cloud ». L'unité de traitement de données 12 peut également être adaptée pour mettre en œuvre des simulations numériques, telles qu'éléments finis, de la pièce à analyser 20. La illustre les étapes du procédé de prévision selon l’invention. Les étapes 800, 801, 802, 804 et 805 utilisent un principe de superposition détaillé ci-dessous. En se référant à la , la pièce fissurée 901 a un FIC K 901 , que l’on cherche à calculer. Selon le principe de superposition formulé par Bueckner, le facteur d'intensité des contraintes de mode I K 901 pour une fissure dans un corps 901 soumis à un champ de contraintes externes 9011 est identique à celui K 903 d'une fissure similaire, soumise à un champ de contraintes internes 9032 dans un corps similaire 903 qui n'est soumis à aucun champ de contrainte externes, où le champ de contrainte interne 9032 agissant sur la fissure est égal et de signe opposé au champ de contraintes 9022 qui existerait à la ligne de fissure dans le corps non fissuré 902 soumis au champ de contraintes externes 9021 égal au champ de contraintes externes 9011 auquel est soumis le corps fissuré 901. Ainsi, on peut réaliser un calcul sur la structure non fissurée 902 et en déduire, par post-traitement du champ de contraintes 9021, les efforts nodaux à appliquer sur les lèvres de fissure de la structure fissurée 903. Afin de calculer le FIC K 901 de la structure fissurée, on utilisera une méthode de calcul connue de l’état de la technique. On peut par exemple utiliser la méthode G-θ ou les intégrales d’interaction (voir référence [7]). Ainsi, les quatre premières étapes du procédé utilisent le principe de Bueckner afin de définir un champ de contraintes internes 9032 à appliquer à la fissure : une étape 800 de maillage de la fissure ; une étape 801 de maillage d’une pièce non fissurée, structurellement identique à la pièce fissurée ; une étape 804 de maillage de la pièce fissurée, par superposition des maillages définis aux étapes 801 et 800 précédentes ; une étape 802 de détermination du champ de contraintes internes 9032. En pratique, le FIC obtenu pour la structure 903 peut être conservatif par rapport au FIC réellement vu par la structure chargée 901, car on impose des conditions aux limites de contraintes 9011 imposées, alors que la structure pourrait par exemple être soumise à des conditions aux limites différentes, notamment des déplacements imposés. Pour de telles conditions, l’avancée de la fissure engendre une perte de raideur de la structure qui fait diminuer les efforts imposés à celle-ci, ainsi que le FIC. En référence à la , et dans le cas où le modèle de calcul comprend un nombre élevé de degrés de liberté et/ou des non-linéarités, par exemple de contact, de comportement, de conditions aux limites variables dues à de grands déplacements dans la structure, ou tout autre non-linéarité, il peut être avantageux d’utiliser des zooms structuraux. On ne conserve alors qu’une partie ou modèle local 1001 du modèle global 1002 dans lequel on souhaite propager la fissure. Ceci permet notamment de réduire considérablement le nombre de nœuds auxquels les calculs devront être effectués, et par conséquent le nombre de valeurs numériques et le temps de calcul nécessaires à l’obtention des FIC ou des durées de vie de la pièce mécanique modélisée. On considère alors les conditions aux limites du modèle local 1001 non-fissuré de sorte qu’elles aboutissent à des champs de contraintes similaires à ceux obtenus aux mêmes points dans le modèle global 1002 non-fissuré. Dans les méthodes de fissuration 3D classiques, cette étape peut être longue, car tant qu’on n’obtient pas des conditions limites du modèle local qui convergent avec les champs de contraintes du modèle global, on doit réitérer le calcul sur le modèle non-fissuré et on ne peut pas commencer le calcul sur le modèle fissuré. En revanche, dans le procédé, l’utilisation du principe de superposition permet de considérer des conditions aux limites simplifiées, à savoir une contrainte imposée sur les lèvres de la fissure, de sorte que le champ de contraintes est connu, et qu’il n’est pas nécessaire de réitérer le calcul. Il peut être avantageux de bloquer les degrés de liberté associés aux mouvements de corps rigides. En effet, de tels mouvements n’interviennent pas dans des calculs de propagation d’une fissure, une telle propagation n’étant affectée que par des efforts internes à la structure et/ou à des variations de température dans la structure. Le procédé selon l’invention présente également l’avantage qu’il n’est pas nécessaire de considérer un champ de température dans la structure, mais qu’on peut spécifier une température uniforme dans la structure. En effet, les effets thermoélastiques qui pourraient avoir un impact sur le calcul des FIC sont contenus dans le champ de contrainte imposé sur les lèvres de la fissure. En règle générale, pour un calcul de propagation de fissure 3D mono-instant, l’instant qui est utilisé dans le calcul, dit instant critique t crit , est celui qui maximise la contrainte principale maximale, c’est-à-dire la contrainte ayant la valeur propre la plus élevée à l’instant de calcul, au nœud représentant le centre de la fissure initiale. Cependant, il est possible que le champ de contrainte ne se comporte pas de manière proportionnelle au passage d’un instant à un autre. Par conséquent, la connaissance du FIC à l’instant critique ne sera pas suffisante pour connaître son évolution aux autres instants par proportionnalité à la contrainte au nœud critique. Une évaluation de la durée de vie en propagation de fissure sur la base de ce seul FIC à l’instant critique pourrait donc se révéler très imprécise. En effet, l’instant critique associé à la maximisation de la contrainte principale maximale au centre de la fissure au cours du temps ne correspond pas nécessairement à l’instant qui maximise le facteur d’intensité des contraintes. L’évolution du FIC en mode I peut être évalué par une équation du type F étant une fonction dépendante de la taille de la fissure (a), de la géométrie de la fissure ( geom ), ainsi que du champ de contrainte sur la surface fissurée normal à la direction principale et normalisé par la contrainte principale au centre de la fissure ( χ ). Le FIC K I dépend donc de cette fonction ainsi que de la contrainte principale au centre de la fissure ( ) et de la taille de la fissure ( a ). Ainsi, de sorte à évaluer la durée de vie d’une pièce de manière précise, il sera nécessaire d’évaluer le FIC sur l’ensemble des instants de la mission effectuée par la pièce. Ceci peut considérablement augmenter le temps de calcul et la quantité de données calculées, rendant une telle méthode rédhibitoire pour des calculs en bureau d’études. De façon à réduire la quantité de données engendrées et le temps de calcul nécessaire à l’obtention des résultats, le procédé selon l’invention comprend une étape 803 de réduction de modèle. Une telle étape permet d’obtenir une évolution temporelle du FIC fiable par pas de propagation de fissure, sans qu’il soit nécessaire de calculer explicitement les FIC à tous les instants et pour chaque pas de propagation de fissure. L’étape de réduction de modèle 803 selon l’invention consiste à modéliser l’évolution temporelle du champ de contrainte par une combinaison linéaire de modes spatiaux. Ces modes sont des champs de contraintes fictifs qui ne dépendent que de l’espace, et qui permettent de reconstituer le champ de contraintes à chaque instant par une combinaison linéaire dépendant exclusivement du temps. Pour chaque pas de propagation simulé, le champ de contraintes ainsi modélisé s’écrit : où n est le nombre de modes nécessaire pour représenter l’évolution temporelle du champ de contrainte avec une précision donnée, σ k est le champ de contraintes correspondant au k-ème mode spatial et a k est la contribution temporelle au k-ème mode spatial, obtenu pour un instant t par le produit scalaire entre le champ de contrainte σ(t i ) et le champ modal σ k . Les x j correspondent aux coordonnées spatiales d’un nœud situé sur les lèvres de la fissure et t i est un instant temporel simulé initialement sur le modèle sain (c’est-à-dire le modèle non fissuré). On obtient ainsi un modèle modal simplifié 805 sur lequel on effectuera une étape 806 de calcul des FIC pour chaque mode structurel, aux lèvres de la fissure. Cette décomposition peut s’obtenir selon plusieurs méthodes : par décomposition orthogonale aux valeurs propres (en anglais Proper Orthogonal Decomposition ou POD), par décomposition en valeurs singulières ( Singular Value Decomposition ou SVD), par décomposition propre généralisée ( Proper Generalized Decomposition ou PGD), par analyse en composantes principales (ACP) ou par toute autre méthode de décomposition appropriée. Le procédé selon l’invention peut ensuite comprendre une étape 807 de reconstruction des valeurs du FIC variant dans le temps à partir des FIC modaux associés aux champs de contrainte spatiaux σ(x j ) , obtenus lors de l’étape 806 de calcul du FIC pour chaque mode structurel : En alternative, il est possible de remplacer les conditions aux limites en contraintes par des conditions aux limites en efforts nodaux. Cette étape 803 de réduction de modèle permet de diminuer considérablement le nombre de calculs à effectuer pour obtenir le FIC : pour un calcul classique du FIC à tous les nœuds du front de fissure et à tous les instants, le nombre de valeurs numériques engendrées par les calculs correspond à l’équation suivante : où T correspond aux nombres d’instants de la mission, I au nombre de nœuds présents dans le maillage de la structure et 6 aux composantes du tenseur des contraintes. En revanche, pour un calcul comprenant l’étape 803 de décomposition, le nombre de valeurs numériques engendrées par les calculs correspond à l’équation suivante : où N correspond au nombre de modes considérés dans le calcul. Ainsi, pour une structure dont le maillage comprend 1000 nœuds, pour laquelle on souhaite évaluer le FIC à 600 instants, et pour lequel on choisit une modélisation selon 5 modes structurels, le nombre de valeurs numériques n c,1 pour une modélisation sans décomposition équivaut à 3 600 000, là où le nombre de valeurs numériques n c,2 pour la modélisation équivalente avec décomposition ne nécessite que 33 000 – soit une réduction de l’espace mémoire nécessaire de plus de 99%. Selon un mode de réalisation du procédé, l’étape de décomposition 803 peut être une décomposition en effort nodal plutôt qu’une décomposition sur le champ de contraintes. L’avantage de raisonner en efforts nodaux plutôt qu’en contraintes est de réduire davantage le nombre de valeurs numériques qui doivent être stockées. En effet, le tenseur des contraintes σ selon l’équation 6 comporte 6 composantes par nœud, alors que le vecteur des efforts F n’en comprend que 3 par nœud. On peut ainsi réduire de moitié le nombre de volume de données numériques nécessaires pour représenter l’intégralité du chargement. Le procédé comprend en outre une étape 808 de détermination du front de fissure a, consistant à calculer l’avancée de la fissure au cours de chaque pas de propagation. Le procédé peut également comprendre une étape 811 de calcul d’une durée de vie d’une pièce mécanique par propagation d’une fissure par fatigue. Les conditions aux limites vues par la structure, notamment de contrainte externe et de température, peuvent évoluer avec le temps. De tels calculs d’avancée de la fissure ainsi que de durée de vie de la pièce mécanique sont effectués dans l’hypothèse que l’évolution de ces conditions aux limites est périodique. Ceci traduit en termes mathématiques le fait que la pièce effectue des missions, notamment des vols dans le domaine aéronautique, et une période correspond ainsi à une mission donnée. Le calcul de durée de vie effectué à l’étape 811 se base sur l’évolution temporelle du FIC obtenue tel que précédemment décrit, ainsi que sur une loi de propagation de fissure. Le procédé selon l’invention pourra se baser sur toute loi de propagation appropriée. Selon un mode de réalisation, l’étape 808 se base sur une loi de Paris. Selon un mode de réalisation, l’étape 808 se base sur une loi de Paris avec correction d’effet de rapport de charge par la loi d’Elber, qui peut s’écrire comme suit : où da/dN est la propagation par fatigue de la fissure en un cycle c de chargement obtenu notamment selon une méthode dite « Rainflow » à laquelle il est fait référence ci-après, C et m sont des coefficients, appelés coefficients de Paris, obtenus de manière expérimentale, ΔKI eff c est l’amplitude du FIC effectif lors d’un cycle, et ΔKI c est l’amplitude de FIC réel lors d’un cycle c , de la méthode dite « Rainflow », T c est la température du cycle c de la méthode dite « Rainflow », R c est le rapport de charge de ce cycle c et a et b sont des coefficients permettant de prendre en compte le fait que la fissure reste fermée en dessous d’un certain FIC. L’utilisation d’une correction d’Elber permet d’obtenir une modélisation de propagation de fissure plus proche de la réalité que l’utilisation d’une simple loi de Paris non corrigée, la correction d’Elber permettant de tenir compte de l’effet du rapport de charge, c’est-à-dire de la fermeture éventuelle de la fissure avant d’atteindre une charge minimale. Selon un autre mode de réalisation, l’étape 808 se base sur une loi de Paris avec correction de Walker. Une telle loi peut s’écrire comme suit : où b w est un facteur représentant l’influence du rapport de charge Rc sur la vitesse de propagation de fissure. L’utilisation d’une correction de Walker permet également d’obtenir une modélisation de propagation de fissure plus proche de la réalité que l’utilisation d’une simple loi de Paris non corrigée, la correction de Walker permettant de tenir compte de l’effet du rapport de charge. Afin d’obtenir l’avancée de fissure à l’étape 808 de détermination du front de fissure de sorte à prendre en compte une multiplicité d’instants, voire tous les instants de la mission modélisée, on pourra notamment utiliser la méthode dite « Rainflow » décrite notamment dans la Norme AFNOR A 03-406. L’avancée de fissure obtenue grâce à l’équation 10 pour chaque point du front de fissure n’est toutefois pas utilisée directement pour la mise à jour de la géométrie de fissure qui initialise le pas de propagation de fissure suivant. En effet, l’avancée de fissure ainsi obtenue est très faible et conduirait à simuler des milliers de pas de propagation numériquement, ce qui est industriellement inenvisageable aussi bien d’un point de vue des temps de calcul que de l’espace mémoire qu’il faudrait allouer pour stocker tous les fichiers de résultats de calcul. Afin de palier à ce problème, il peut être prévu qu’on spécifie une valeur d’avancée de fissure maximale a max que l’on impose pour faire propager la fissure. Cette valeur a max varie avec la propagation de la fissure, on aura donc une valeur a max distincte pour chaque instant de la modélisation. On peut ainsi réduire le nombre de pas de propagation à simuler. En contrepartie on ne dispose plus des informations concernant les pas de propagation intermédiaires. L’équation 10 caractérisant la loi de propagation sera par exemple ainsi modifiée : où la valeur du coefficient est fixée de telle sorte que le point du front de fissure qui se propage le plus par la loi selon l’équation 10 se retrouve avec une avancée de fissure maximale (au point du front de fissure qui avance le plus) égale à amax. Cette opération constitue un dé-raffinement de la discrétisation du calcul en pas de propagation, et permet ainsi une réduction du nombre de pas de calcul de l’étape 808 d’estimation de la valeur d’avancée de fissure. Selon un mode de réalisation, le procédé peut prendre en compte une contribution à l’avancée de fissure due au temps de maintien. En effet, l’avancée de fissure peut être plus importante pour un cycle donné, lorsque la pièce mécanique est maintenue à un FIC et/ou à une température suffisamment élevée pendant un certain temps, quand bien même la pièce n’est pas soumise à une variation du FIC généralement associée au phénomène de fatigue. La propagation de la fissure s’arrête sous deux conditions. Premièrement, lorsque la fissure a traversé toute la pièce. Deuxièmement, lorsque le FIC mode I, K I , atteint la ténacité du matériau, la ténacité étant une mesure de la résistance du matériau considéré à la propagation de fissures indépendamment de la structure considérée. Comme mentionné précédemment, l’étape de détermination 808 du front de fissure se base sur une direction de propagation imposée, qui est généralement celle qui maximise la contrainte principale. L’angle de bifurcation par rapport au plan initial de la fissure, et pour chaque nœud du front de fissure, est alors donné par l’équation suivante (on utilise les facteurs d’intensité des contraintes avant propagation) : où les FIC K I , K II utilisés sont ceux avant propagation de la fissure. Lors de l’étape de reconstitution (807) du facteur d’intensité des contraintes (FIC) à tous les instants de la modélisation, les facteurs d’intensité des contraintes associés aux modes I, II et III peuvent varier temporellement sans corrélation aucune entre eux, c’est-à-dire que l’instant qui maximise chaque FIC peut varier d’un nœud à l’autre. Ainsi, une approximation qui consisterait à prendre en compte, pour chaque nœud, l’instant qui maximise le FIC à ce nœud mènerait à des résultats trop conservatifs. Il est donc nécessaire, pour appliquer un tel critère, de déterminer pour chaque nœud, à quel instant il faut récupérer les facteurs d’intensité des contraintes associés aux modes I et II ( K I et K II ) pour pouvoir utiliser l’équation 12. Selon un mode de réalisation, le procédé selon l’invention utilise à cet effet des formules issues de travaux théoriques portant sur la mécanique de la rupture, de sorte à obtenir les trois FIC K I , K II et K III . L’état de la technique enseigne qu’il est possible de déduire l’évolution du facteur d’intensité des contraintes après propagation en fonction d’une petite extension de fissure δ (valeur maximale de l’avancée du front de fissure) et ce, même si la fissure a bifurqué. On peut notamment utiliser l’équation : où les termes à gauche de l’équation sont les FIC après extension de la fissure, et les termes KI*, KII* et KIII* sont les FIC immédiatement après bifurcation de la fissure. Tous les autres termes de la somme (à savoir ), sont négligeables pour une extension de fissure δ suffisamment petite. Le terme δ.η(x 3 ) donne la forme de l’extension, et la valeur de l’extension. x 3 correspond à l’abscisse curviligne le long du front de fissure. Selon un mode de réalisation préférentiel, le procédé selon l’invention propose une autre méthode de gestion de l’orientation de la fissure. Cette méthode consiste à rechercher l’angle qui maximise l’avancée de fissure da/dN lorsque l’on remplace le FIC de mode I avant propagation ( K I ) par le FIC de mode I après propagation ( K I * ) dans la loi de propagation, par exemple définie selon l’équation 12. Ceci permet notamment de mettre en œuvre le calcul de l’avancée de fissure détaillé précédemment pour chaque point du front de fissure et pour différents angles de bifurcation. Pour chaque nœud du front de fissure donné, il suffit ensuite de chercher l’angle de bifurcation φ, parmi une série de valeurs représentatives d’une plage choisie, qui maximise da/dN . On note cet angle de bifurcation φ a et la vitesse de propagation correspondante ( da/dN) a . Cet angle n’est pas forcément optimum dans la mesure où une incertitude peut lui être associée. Les valeurs de φ a et de ( da/dN ) a vont servir de valeurs initiales pour un algorithme d’optimisation. L’algorithme d’optimisation se base ainsi sur ces deux valeurs, ainsi que sur leurs équivalents aux angles précédant et suivant dans la suite d’angles de bifurcation recherchés, notés respectivement φ a-1 et φ a+1 ainsi que les vitesses de propagation associées ( da/dN ) a-1 et ( da/dN ) a+1 Ces différentes valeurs vont servir à identifier un polynôme de degré 2 ( α , ß et γ étant les coefficients du polynôme) par la méthode dite des « moindres carrés » : Chercher la valeur analytique de φ pour laquelle la dérivée de P( φ) par rapport à φ s’annule permet de trouver une valeur potentielle de φ maximisant la vitesse de propagation, de façon plus précise que φ a . Notons cette valeur φ a * et la vitesse de propagation associée (da/dN) a * obtenue par une combinaison des formules selon les équations 11 et 14. On dispose désormais de 4 angles et de 4 vitesses de propagation. On ré-identifie le polynôme selon l’équation 15 avec les angles associés aux trois vitesses de propagation les plus élevées et on recommence jusqu’à aboutir à une convergence. Si la valeur de φ qui annule le polynôme varie de moins d’une valeur prédéterminée, par exemple mais non limitativement 1/1000, par rapport à la valeur de φ a * de l’itération précédente, alors l’algorithme s’arrête. On considèrera que la valeur de φ a * ainsi obtenue est optimale pour le nœud du front de fissure considéré. La méthode de gestion de l’orientation de la fissure ainsi définie n’utilise que le facteur d’intensité des contraintes associé au mode I. Cependant, le procédé pourra utiliser une méthode étendue à la prise en compte des modes II et III par le biais d’un facteur d’intensité des contraintes équivalent s’exprimant comme une combinaison de ceux en modes I, II et III. La méthode de gestion de l’orientation du procédé selon l’invention présente l’avantage de ne comporter que des calculs analytiques, contrairement aux méthodes de fissuration 3D selon l’état de la technique, qui nécessitent par exemple des calculs d’éléments finis. Ceci permet notamment de réaliser les calculs pour chaque point du front, pour chaque instant de la mission et à chaque pas de propagation, sans impact majeur sur le coût et le temps de calcul de la méthode. On peut considérer que les pas de propagation obtenus par le procédé sont suffisamment faibles pour pouvoir simplifier l’équation 14 en : où φ(x 3 ) correspond à l’angle de bifurcation. Les termes de la matrice de couplage, à savoir les termes F 1,1 ; F 1,2 ; F 1,3 ; F 2,1 ; F 2,2 ; F 2,3 ; F 3,1 ; F 3,2 et F 3,3 sont donnés dans la littérature : F i,j ( φ (x 3 )), avec i et j étant 1 ou 2 sont donnés dans les références (1) et (2) ; F i,3 ( φ (x 3 )) et F 3,i ( φ (x 3 )), avec i étant 1 ou 2, sont donnés dans les références (1), (2) et (3) ; F 3,3 ( φ (x 3 )) est donné dans les références (1) et (3). Il est ainsi possible d’exprimer les FIC après extension de fissure et après bifurcation avec suffisamment de précision tant que l’extension de fissure n’est pas trop élevée, ce qui est généralement le cas. Selon le procédé, on peut déterminer l’angle de bifurcation φ pour un nœud du front de fissure donné à l’abscisse curviligne x 3 , en cherchant l’instant et l’angle qui maximisent K I * (x 3 ) . Ceci permet de généraliser l’équation 13 dans les cas où la direction de propagation diffère pour chaque nœud du front. On définit ainsi une approche plus conforme aux résultats d’essais. En effet, si le chargement consiste en un nombre de cycles conséquent, il se peut que le cycle principal (pour lequel la valeur maximale du facteur d’intensité des contraintes correspond à la contrainte principale maximale à l’instant critique t crit au nœud du centre de la fissure) conduise à une avancée de fissure inférieure à celle induite par le cumul des cycles secondaires associés à des champs de contraintes qui ne sont pas proportionnels à celui associé à la contrainte principale maximale à l’instant critique au nœud du centre de la fissure. Pour un pareil cas, on s’attend à ce que l’angle de bifurcation soit plutôt majoritairement dépendant du chargement associé aux cycles secondaires. Le critère de bifurcation à l’instant critique, donné à l’équation 13, n’est alors pas du tout applicable. À l’issue des étapes de reconstitution du facteur d’intensité des contraintes à tous les instants de la modélisation 807 et de détermination du front de fissure 808, on dispose de valeurs des FIC K I , K II et K III pour chaque nœud du front de fissure, ainsi que des coordonnées spatiales des fronts de fissure. On peut ainsi reconstituer, selon une étape d’obtention d’un maillage fissuré propagé 809, un maillage de la structure avec la fissure propagée. Le procédé peut ensuite comprendre une étape 810 d’établissement de tableaux donnant pour chaque dimension de fissure et pour chaque instant, une valeur de facteur d’intensité des contraintes. Ces tableaux sont appelés « formulaires ». Enfin, les tableaux établis à l’étape 810 d’établissement de formulaires seront utilisés comme données d’entrée dans un solveur de calcul de la durée de vie en propagation de fissure selon une étape 811 de calcul de la durée de vie de la pièce mécanique. L’étape 810 d’établissement de formulaires ne prend en compte que les valeurs du FIC en mode I, K I. En effet, les codes de calcul de durée de vie utilisés à l’étape 811 de calcul de la durée de vie prennent en entrée une information unidirectionnelle. On post-traite alors les résultats issus de la propagation de fissure 3D pour se rapporter à une information unidirectionnelle. À cet effet, l’étape 810 d’établissement de formulaires comprend la réduction dimensionnelle d’un calcul 3D complexe à un simple tableau 1D. Selon un mode de réalisation, cette réduction dimensionnelle de l’étape 810 d’établissement de formulaires se fait en ramenant la problématique de la géométrie de la fissure en trois dimensions à une problématique de fissure plane à front droit, qui est un modèle standard bien connu de l’homme du métier. Cette réduction, représentée en et b, consiste alors à évaluer une équivalence thermodynamique fondée sur une égalité de l’énergie dissipée entre les deux modélisations. Le fondement physique est expliqué ci-dessous (voir aussi référence [5]). Le premier principe de la thermodynamique, pour une transformation du front de fissure infinitésimale donne, exprimée en densité d'énergie par unité de longueur de front de fissure : Le second principe de la thermodynamique permet d'exprimer l'irréversibilité de la transformation : La densité d'énergie dissipée est toujours positive. Le terme da représente l'incrément de surface fissurée par unité de longueur de front de fissure au cours de la transformation infinitésimale. Cette quantité est homogène à une longueur : il s'agit en quelque sorte d'un incrément de longueur de fissure au sens « thermodynamique » du terme. Il apparaît ainsi une grandeur « physique » que l'on pourra considérer pour effectuer une interpolation qui présente un sens physique. Les éléments précédents permettent de rendre le problème de fissuration tridimensionnel initial équivalent à un problème de propagation de fissure plane à front droit en élasticité plane. Pour plus de généralité, nous nous plaçons dans le cas de l'élasticité adaptée, c'est-à-dire que le matériau a déjà pu subir une déformation plastique initialement, mais la propagation de fissure se fait sous un chargement cyclique de fatigue en élasticité. Le champ de température peut évoluer au cours du temps et est supposé ne pas varier spatialement. On considère en général, la température au nœud du maillage initial (non-fissuré) correspondant au centre de la fissure initial injectée dans le modèle. Par conséquent, à un cycle de fatigue est associée une température. Ainsi, la réduction dimensionnelle se fait de manière à respecter la conservation de certaines contraintes physiques : Conservation des énergies dissipées par cycles et entre chaque pas de propagation de fissure issu du calcul de fissuration 3D ; Alternativement au point précédent, conservation de la longueur des fronts de fissures et des surfaces fissurées ; Conservation de la cohérence en termes de propagation de fissure vis-à-vis de la loi de propagation de fissure utilisée. On pourra utiliser une loi de propagation telle que celle formulée à l’équation 12. Alternativement, on pourra utiliser une autre loi de propagation, telle qu’une loi de Paris non corrigée ou une loi plus complexe. La mise en œuvre de cette méthode de réduction dimensionnelle est du type de celle décrite dans la demande de brevet WO2017/220923, à laquelle on pourra avantageusement se référer, et conduit à l’établissement de deux formulaires donnant le FIC KI équivalent en fonction de la taille de fissure. L’un est obtenu avec la méthode dite « énergétique » ( ) et l’autre avec la méthode dite « géométrique » ( ). En référence à la , qui illustre la méthode énergétique, le haut de la figure représente le modèle tridimensionnel, où le front de fissure se propage, sur deux pas de propagation, d’un front 23 à un front 22 puis à un front 21. Une énergie 24 est dissipée durant la propagation de la fissure du front 23 au front 22, et une énergie 25 est dissipée durant la propagation de la fissure du front 22 au front 21. On passe à un modèle unidimensionnel en conservant l’énergie dissipée, de telle sorte que la dissipation unidimensionnelle 27 sera égale à la dissipation tridimensionnelle 24, et que la dissipation unidimensionnelle 28 sera égale à la dissipation tridimensionnelle 25. En référence à la , qui illustre la méthode géométrique, le haut de la figure représente le modèle tridimensionnel, où le front de fissure se propage, sur deux pas de propagation, d’un front 23’ à un front 22’ puis à un front 21’. Une surface 24’ de la pièce est fissurée durant la propagation de la fissure du front 23’ au front 22’, et une surface 25’ est fissurée durant la propagation de la fissure du front 22’ au front 21’. On passe à un modèle unidimensionnel en conservant la surface fissurée, de telle sorte que la surface fissurée unidimensionnelle 27’ sera égale à la surface fissurée tridimensionnelle 24’, et que la surface fissurée unidimensionnelle 28’ sera égale à la surface fissurée tridimensionnelle 25’. Ces deux méthodes de réduction dimensionnelle peuvent être utilisées dans un code de calcul de la durée de vie en propagation de fissure (tel que par exemple le code AFGROW développé par les laboratoires de recherche de l’armée de l’air des Etats-Unis). La durée de vie retenue sera la plus faible des deux durées de vie obtenues par ces deux formulaires. Dans le cas où le calcul de fissuration 3D est mono-instant, c’est-à-dire calculé à l’instant critique, avant l’étape 807 de reconstitution du facteur d’intensité à tous les instants, ce formulaire est construit au niveau de l’instant critique. On reconstruit alors les formulaires pour les autres instants par proportionnalité entre le K I équivalent et la contrainte principale maximale au niveau de l’anomalie. La taille de fissure équivalente discrétisée ne change pas d’un instant à un autre. Dans le cas où le calcul de fissuration 3D est multi-instants, c’est-à-dire qu’on se trouve après l’étape 807 de reconstitution du facteur d’intensité de contraintes à tous les instants, on dispose de valeurs de FIC pour chaque point des fronts de fissures et pour chaque instant de la mission. Il est alors possible de construire des formulaires pour chaque instant de la mission pour les fournir aux logiciels de calcul de durée de vie en utilisant les équations 6 et 9. Selon un autre mode de réalisation, la réduction dimensionnelle de l’étape 810 d’établissement des formulaires se fait en considérant des directions partant du centre de la fissure initiale et en les projetant sur les fronts de fissures. Pour chaque direction on récupère les distances entre les fronts successifs au niveau des points de projection sur les fronts avec la direction donnée ainsi que les facteurs d’intensité des contraintes sur les points de projection. On détermine la taille de fissure pour un front de fissure donné comme la somme des distances inter-fronts en partant du centre de la fissure initiale. Ainsi, on construit pour chaque direction un formulaire donnant le facteur d’intensité des contraintes en fonction de la dimension de fissure. La durée de vie retenue sera la durée de vie minimale obtenue avec les formulaires. Comme pour le mode de réalisation précédemment décrit, la réduction dimensionnelle peut se faire pour le calcul mono-instant ou pour le calcul multi-instants. Dans le cas multi-instants, on dispose de valeurs de FIC pour chaque point des fronts de fissures et pour chaque instant de la mission. Il est alors possible de construire des formulaires pour chaque instant de la mission pour les fournir aux logiciels de calcul de durée de vie. Une fois les formulaires établis pour un nombre de longueurs de fissure données, pour le calcul mono-instant ou pluri-instants, l’étape 810 d’établissement de formulaires peut comprendre l’établissement de formulaires comprenant un plus grand nombre de longueurs de fissures avec un pas moins élevé par interpolation selon une loi racine. En d’autres termes, on construit un autre formulaire fondé sur le premier formulaire, donnant le carré du FIC en fonction de la dimension de la fissure. A partir de ce tableau, il est possible d’ajouter des longueurs de fissures intermédiaires par interpolation linéaire des carrés des FIC mode I K I présents dans les formulaires non interpolés. Il est ensuite possible de reconstruire le tableau K I vs taille de fissure interpolée. La décision d’ajouter des points se fait par comparaison des valeurs de K I ajoutées obtenues par les méthodes d’interpolation linéaire sur le K I d’une part et sur le carré du K I d’autre part comme il sera expliqué plus tard. Si les écarts entre ces deux valeurs dépassent un seuil prédéfini, par exemple mais de manière non limitative 0,3%, le point est ajouté. On ajoute autant de points qu’il faut pour qu’il ne puisse pas y avoir d’écart excédant le seuil prédéfini sur toute la plage de valeur décrite par le formulaire interpolé donnant le facteur d’intensité des contraintes en fonction de la taille de fissure. Cette interpolation effectuée sur les carrés des K I permet d’optimiser l’espacement des points du formulaire, de sorte à améliorer la précision numérique sur le calcul de la durée de vie en propagation de fissure. Avantageusement, et afin de limiter les coûts de calcul, le procédé ne comprend pas le calcul de l’avancée et de la direction de fissure induites par une mission de chargement complète, mais détermine plutôt une valeur d’avancée maximale de fissure a max initiale, c’est-à-dire pour une première propagation, de sorte que pour une fissure de taille initiale a 0 , le passage de a 0 à a 0 + a max est convergé en a max . Il est avantageux de choisir un a max optimal, ni trop petit pour éviter des coûts de calcul trop élevés, ni trop grand pour éviter des calculs imprécis correspondant à une perte d’information. L’évolution du a max au cours de la propagation peut être automatisée en appliquant les étapes suivantes : Considérer la valeur de a max constante pour un nombre P de premiers pas, par exemple mais non limitativement pour 5 pas. Pour les pas suivants, une optimisation de la valeur de a max est possible. On procède comme suit. On établit deux interpolations du carré des facteurs d’intensité des contraintes effectifs ΔK eff 2 en fonction de la dimension de la fissure équivalente. La première est une interpolation fondée sur une courbe de régularité C 2 . La deuxième consiste en une interpolation linéaire des points. On appelle δ l’écart maximal entre ces deux interpolations sur l’intervalle considéré. Si cet écart δ est inférieur à un certain seuil prédéterminé, par exemple mais non limitativement 0.3%, pour un nombre d’intervalles consécutifs précédents, alors la valeur de a max est augmentée pour le pas de propagation suivant, par exemple mais non limitativement de 10% A l’inverse, si l’écart δ est supérieur au seuil prédéterminé sur l’intervalle précédent, on réalise un pas de propagation 3D intermédiaire puis la valeur de a max est réduite, pour le pas de propagation suivant, par exemple mais non limitativement de 10%. Ainsi, à l’issue de l’étape 810 d’établissement des formulaires, on obtient un formulaire final donnant le facteur d’intensité des contraintes en fonction de la dimension de la fissure. L’étape 810 d’établissement des formulaires peut s’appliquer au calcul mono-instant, ou au calcul multi-instants à chaque instant modélisé, postérieurement à l’étape 807 de reconstitution du facteur d’intensité des contraintes à tous les instants. En référence à la , et afin de réduire davantage le nombre de valeurs numériques à obtenir et le temps nécessaire à la mise en œuvre du procédé, il peut être avantageux de concentrer l’effort de résolution du problème numérique, résolu par exemple par la méthode des éléments finis, sur une zone remaillée 32, qui n’est qu’une partie de la zone globale modélisée 30. À cet effet, le procédé comprend, selon un mode de réalisation, une étape additionnelle de décomposition de la structure en deux zones. On peut décrire le problème numérique à résoudre de la forme suivante : où s’interprète comme étant la matrice de rigidité, étant les valeurs nodales associées au champ de déplacement et étant le vecteur second membre contenant les informations nodales sur les conditions aux limites. Les matrices , et peuvent se décomposer sous la forme suivante : Où l’indice r fait référence aux degrés de liberté retenus et l’indice n fait référence aux degrés de liberté non retenus. Avec cette notation, l’équation 19 devient : La première du système d’équations précédent conduit à : Ceci permet de ré-exprimer les degrés de liberté de à partir de ceux de qui nous intéressent davantage. Avec la deuxième équation du système d’équations 21, on a : D’où l’on déduit : Et donc on peut se ramener à un système matriciel très similaire au système linéaire initial, mais sur un nombre de degrés de liberté plus faible : Cette technique peut s’appliquer en fissuration 3D. En effet, l’insertion d’une fissure dans un maillage initial qui n’en contenait pas ne conduit pas, généralement, à un remaillage intégral de la structure. Seuls les degrés de liberté de la zone remaillée 32 interviennent dans le calcul du facteur d’intensité des contraintes. Ainsi, les degrés de liberté d’indice r peuvent être les degrés de libertés associés à la zone remaillée 32, tandis que ceux d’indice n peuvent être ceux de la zone non-remaillée 31. Ceci permet avantageusement de résoudre un problème numérique (équation 25) ne comprenant que les degrés de liberté correspondant à la zone remaillée 32, plutôt qu’un problème numérique (équation 19) correspondant à l’intégralité des degrés de liberté de la structure modélisée 30. La zone non-remaillée 31 le restera pour plusieurs pas de propagation. Ainsi, tant que la zone non-remaillée ne change pas, peut être calculée et inversée ( ) une seule fois, par exemple au moyen de la méthode de Cholesky. Le calcul de peut être plus ou moins long selon que l’on ait recours ou non à un zoom structural tel que mentionné précédemment. Avantageusement, le calcul de ne se fait qu’une seule fois au cours de la simulation de fissuration 3D. Ceci permet de réduire le nombre de calculs, surtout lorsque le nombre de pas de propagation à simuler est élevé. Si toutefois on souhaite effectuer le calcul pour les degrés de liberté d’indice n, on pourra utiliser l’équation 21. L’équation 24 peut être résolue par n’importe quel solveur du commerce, par exemple mais non limitativement via une décomposition LU ou via une méthode itérative. Une telle décomposition de la structure en une zone remaillée 32 et une zone non-remaillée 31 permet avantageusement de réduire les temps de calcul, notamment pour pouvoir considérer l’ensemble de la structure fissurée plutôt qu’un zoom structural tel que défini précédemment, qui pourrait artificiellement assouplir la structure. Cette étape permet donc avantageusement d’obtenir un calcul de durée de vie moins conservatif, plus proche de la durée de vie réelle de la structure, et permet donc de se rapprocher voire de satisfaire les objectifs attendus de durée de vie pour une pièce mécanique. Il convient de noter que, dans le cas d’un calcul multi-instants, les matrices et ne dépendent que des propriétés matériaux. On peut supposer qu’elles restent inchangées pour les différents modes. Elles sont donc indépendantes des modes issus de décomposition opérée à l’étape 803 de réduction de modèle. Il suffit donc de n’en stocker qu’une seule copie, correspondant à un des différents modes. En revanche, les matrices et dépendent explicitement des modes issus de la décomposition opérée à l’étape 803 de réduction de modèle. Il faudra donc les calculer pour chaque mode. Par conséquent, les vecteurs et dépendent également de ces modes. C’est à partir des et des que les FIC modaux pourront être calculés. Dans toute la description, le procédé selon l’invention comprend un calcul des facteurs d’intensité de contraintes K I . Cependant, pour des matériaux non-homogènes et/ou non isotropes, la notion de FIC est mal définie, de sorte qu’on pourra remplacer les FIC par leur taux de restitution d’énergie G I équivalent. Le procédé selon l’invention peut présenter des applications avantageuses, notamment dans le domaine de l’aéronautique. La durée de vie d’une pièce mécanique estimée grâce au procédé de prévision selon l’invention peut permettre, une fois l’apparition d’une anomalie - modélisée par une fissure de géométrie similaire-, constatée dans une pièce donnée, de calculer combien de missions en vol cette pièce pourra effectuer dans le futur. Alternativement, on pourra utiliser le procédé de prévision de manière à calculer, pour différentes géométries et/ou positionnements initiaux de fissure dans la pièce, combien de missions en vol cette pièce pourra effectuer. En référence à la , l’invention porte également sur un procédé de maintenance d’une pièce d’un aéronef entre deux séries M n et M n+1 de missions en vol. Chaque série peut comprendre une unique mission en vol ou une pluralité de missions, on définira alors les valeurs ΔM n et ΔM n+1 comme le nombre de missions de la première et de la deuxième série de missions respectives. La première série de missions M n peut éventuellement être précédée d’autres missions en vol, et la deuxième série de missions M n+1 peut quant à elle être suivie d’autres missions en vol. Le procédé de maintenance est effectué sur un aéronef au sol, entre la première et la deuxième série de missions. Un opérateur ou une machine mettra en œuvre une étape d’inspection 110 de la pièce, par exemple mais de manière non limitative d’inspection visuelle, afin de déterminer si une ou des anomalie(s) sont apparues dans la pièce à inspecter suite à la première série M n de missions en vol. Si aucune anomalie n’est détectée, une étape de remise en service 112 de l’aéronef pourra être effectuée sans remplacement ou réparation de la pièce. Si une anomalie est détectée, on mettra en œuvre le procédé de simulation défini précédemment afin de calculer une durée de vie restante de la pièce, en modélisant l’anomalie par une fissure de géométrie similaire. Si cette durée de vie est supérieure au nombre de missions ΔM n+1 restant dans la deuxième série M n+1 de missions, alors on pourra également mettre en œuvre une étape de remise en service 112 de l’aéronef, sans remplacement ou réparation de la pièce. En revanche, si la durée de vie obtenue est inférieure ou égale au nombre de missions ΔM n+1 restant dans la deuxième série M n+1 de missions, alors on effectuera une étape 111 de réparation ou de remplacement de la pièce, ou de simple mise hors service de l’aéronef selon les cas. L’invention porte en outre sur un procédé d’établissement d’une base de données concernant une pièce mécanique. Ce procédé est mis en œuvre par un système comprenant une unité de traitement de données 10. On pourra ainsi établir une base de données mettant en relation des données d’entrée relatives à la pièce, à savoir des champs de contraintes externes, des tailles et/ou géométries de fissure et des positionnements de fissure dans la pièce, ainsi que les FICs et/ou les durées de vie obtenus pour ces données d’entrée au moyen du procédé de prévision précédemment détaillé. Cela permettra notamment, avant même l’entrée en service d’une pièce mécanique, de pouvoir disposer d’informations précises sur l’évolution d’une fissure et sur la durée de vie de la pièce selon différents scénarios. Procédé de prédiction de la propagation d’une fissure par fatigue dans une pièce mécanique soumise à un champ de contraintes externe (9011) et modélisée numériquement entre un instant initial t 0 et un instant final t f avec un pas de Δt, comprenant les étapes suivantes : une étape de maillage à l’instant t 0 ; une étape de détermination (802) d’un champ de contraintes internes (9032) à appliquer aux lèvres de la fissure, par calcul du champ de contraintes (9022) aux mêmes points du maillage pour la structure non fissurée (902) lorsque le champ de contraintes externes (9011) lui est appliqué, le champ de contrainte interne (9032) de la pièce réelle fissurée (903) étant alors égal et de signe opposé au champ de contrainte interne (9022) de la pièce non fissurée équivalente (902) ; caractérisé en ce que l’étape de maillage comprend les sous-étapes suivantes : Maillage (801) d’une pièce non fissurée (902) identique à la pièce fissurée (901) ; Maillage numérique (800) de la fissure ; Insertion du maillage de fissure (800) dans le maillage (801) de la pièce non fissurée (902), de sorte à obtenir un maillage numérique (804) de la pièce fissurée (901) ; Le procédé comprenant en outre une étape de réduction de modèle (803) de la structure réelle à un modèle modal simplifié (805), la réduction de modèle consistant à exprimer le champ de contrainte à chaque nœud du maillage fissuré (804) sous forme d’une combinaison linéaire de champs de contraintes modaux (σ k ) associés chacun respectivement à un mode structurel k de la pièce, pour un nombre de modes prédéterminé ; une étape de calcul de facteurs d’intensité des contraintes modaux (806) le long d’un front de fissure pour le modèle modal simplifié (805) ; une étape de reconstitution (807) de l’évolution temporelle des FIC à chaque nœud du front de la fissure à partir des FIC modaux obtenus à l’étape (806) de calcul ; une étape d’estimation (808) de la valeur d’une avancée de fissure à partir des évolutions temporelles des FIC obtenues à l’étape (807) de reconstitution et au moyen d’une loi de propagation de fissure, l’étape d’estimation (808) comprenant la détermination de la direction d’avancée de chaque nœud du front de fissure par calcul de l’angle qui maximise la vitesse de propagation de ce nœud ; une étape (809) d’obtention d’un maillage fissuré propagé reconstitué à partir du maillage (804) de la pièce fissurée et de l’estimation (808) de la valeur d’avancée de fissure ; les étapes b. à g. étant répétées pour chaque pas de la propagation modélisée, de l’instant t 0 + Δ t jusqu’à l’instant final t f , le procédé étant mis en œuvre par une unité de traitement de données (12) d’un système (10), l’unité de traitement comprenant des moyens de calcul (14). Procédé selon la revendication précédente, comprenant en outre : une étape (810) d’établissement de formulaires consistant à établir, pour chaque pas de propagation, des tableaux de correspondance entre des facteurs d’intensité de contraintes calculés et des dimensions de fissure ; une étape (811) de calcul d’une durée de vie restante de la pièce mécanique sur la base des formulaires établis. Procédé selon l’une quelconque des revendications précédentes, la réduction de modèle (803) étant une décomposition orthogonale propre (POD). Procédé selon l’une quelconque des revendications précédentes, la loi de propagation de fissure utilisée dans l’estimation (808) de la valeur d’avancée de fissure étant une loi de Paris. Procédé selon la revendication précédente, la loi de propagation de fissure utilisée dans l’estimation (808) de la valeur d’avancée de fissure étant une loi de Paris à laquelle une correction est appliquée. Procédé selon l’une quelconque des revendications précédentes, l’étape d’estimation (808) de la valeur d’avancée de fissure comprenant l’imposition, à chaque pas de propagation, d’une valeur de propagation maximale a max de la fissure de façon à réduire le nombre de pas de propagation calculés pour cette étape. Procédé selon l’une quelconque des revendications précédentes, le maillage de la structure fissurée étant décomposé en : une zone remaillée (32), à proximité de la fissure, pour laquelle les étapes d’obtention d’un maillage numérique (800) de la fissure et d’obtention d’un maillage numérique (804) de la pièce fissurée sont renouvelées à chaque pas ; une zone non-remaillée (31), en dehors de la zone remaillée (32), pour laquelle le maillage est déterminé à l’étape d’obtention d’un maillage numérique (801) de la pièce non fissurée et est maintenu pour l’ensemble des pas de propagation de la fissure ; et en ce que les étapes de calcul de facteurs d’intensité de contrainte (806) le long du front de fissure pour chaque mode structurel du modèle simplifié, d’obtention (809) d’un maillage fissuré propagé et le cas échéant de reconstitution (807) du facteur d’intensité des contraintes à tous les instants de la modélisation ne sont effectuées que pour les degrés de liberté de la pièce correspondant à la zone remaillée (32). Procédé de maintenance d’une pièce d’aéronef entre une première série (M n ) de missions en vol et une deuxième série (M n+1 ) de missions en vol, comprenant les étapes suivantes : Inspection (110) de la pièce afin de déterminer si une anomalie susceptible de causer l’apparition d’une fissure est apparue dans la pièce suite à la première série (M n ) de missions en vol Si oui, mise en œuvre du procédé selon l’une quelconque des revendications 1 à 7 afin de déterminer une durée de vie de la pièce ; Sinon, remise en service (112) de l’aéronef pour la deuxième série de missions en vol ; Procédé de maintenance selon la revendication précédente, dans lequel l’étape a.i. de détermination d’une durée de vie est suivie d’une comparaison avec une durée de vie de la pièce non fissurée et/ou de comparaison avec une durée de vie minimum nécessaire à une certification de la pièce. Procédé d’établissement d’une base de données pour une pièce mécanique, reliant des données d’entrée à des données de sortie ; les données d’entrée comprenant des tailles et/ou géométries de fissures, des positionnements de fissures dans la pièce et des contraintes externes à la pièce ; les données de sortie comprenant des facteurs d’intensité des contraintes et/ou des durées de vie de la pièce, comprenant les étapes suivantes : Mise en œuvre du procédé selon l’une quelconque des revendications 1 à 9 de sorte à obtenir des facteurs d’intensité des contraintes et/ou une durée de vie de la pièce mécanique pour les conditions d’entrée ; Enregistrement, dans la base de données, d’une correspondance entre les données d’entrée et les données de sortie. Système comprenant une unité de traitement (12) comprenant des moyens de calcul (14) et une mémoire (16), l’unité étant configurée pour mettre en œuvre un procédé selon l’une quelconque des revendications précédentes. Produit programme d’ordinateur configuré pour être mis en œuvre par un système selon la revendication précédente, et comprenant des instructions pour provoquer la mise en œuvre d’un procédé selon l’une quelconque des revendications 1 à 10.