Une méthode d'éléments finis adaptative pour les écoulements instationnaires complexes

Nous présentons dans cette thèse le développement d'une procédure générale de résolution des écoulements instationnaires de fluides Newtoniens, construite à partir de la méthode des éléments finis. Ces travaux s'identifiant toutefois davantage à une preuve de concept, la procédure a été formulée dans un contexte 2D uniquement. Cette procédure a pour objectif général de répondre à la difficulté rarement mentionnée de simuler tant adéquatement qu'efficacement cette classe d'écoulements. Pour ce faire, des méthodes adaptatives en espace (adaptation de maillage) et en temps (intégration temporelle adaptative), dont la pertinence a été longuement éprouvée au cours de travaux de recherche antérieurs, sont mises à profit. Se faisant, nous entendons apporter des garanties de précision sur les résultats obtenus, et optimiser le coût de calcul en cours de simulation. La procédure adaptative ainsi mise en place est complètement automatique. Aucune intervention de l'utilisateur n'est donc requise entre l'instant initial et final de l'intervalle de temps simulé. Des procédures similaires ont déjà été proposées dans la littérature, mais l'originalité de celle exposée dans ces pages réside dans le choix et le fonctionnement de ses composantes. Nous avons également mis un accent tout particulier sur la question du transfert d'une solution éléments finis entre deux maillages. Cette opération intervenant après chaque adaptation de maillage, il convient de pouvoir contrôler les inévitables erreurs qui en résultent. L'adaptation du maillage, combinant une méthode de type h et une méthode de type r, permet d'ajuster localement la finesse du maillage au cours de la simulation, tout en préservant sa qualité d'ensemble. La méthode de type h repose sur une estimation de l'erreur spatiale (commise sur chaque élément du domaine de calcul), déduite d'une reconstruction des gradients de la solution éléments finis sur des patchs nodales. Selon le nombre de variables dépendantes actives pour le problème simulé, l'erreur spatiale peut s'exprimer sous la forme d'une ou plusieurs norme(s) d'erreur (approche mononorme et multinorme, respectivement). L'usage d'un opérateur de transition, dont la formulation vise une équidistribution de l'erreur élémentaire pour chaque norme d'erreur, permet de déterminer la taille optimale de chaque élément du maillage. Qui plus est, nous avons utilisé la formulation de l'opérateur de transition pour identifier les éléments du maillage les plus problématiques, soient ceux dont la taille optimale est sensiblement différente de leur taille courante. En retirant uniquement du maillage courant ces éléments et leurs voisins, nous générons ainsi des poches (ou cavités) que nous discrétisons à nouveau avec des éléments d'une taille plus adaptée. Un mailleur frontal est utilisé pour conduire cette dernière opération. La méthode d'adaptation de maillage de type r, résultant de la combinaison de plusieurs méthodes de relaxation de maillage, est finalement appliquée pour maintenir la qualité d'ensemble des éléments du maillage. Ici encore, nous proposons une formulation locale de cette deuxième méthode d'adaptation, en nous appuyant sur une métrique de qualité élémentaire bien choisie. La relaxation du maillage peut ainsi être réalisée uniquement dans les régions du domaine qui le nécessitent. Dans l'ensemble nous proposons et utilisons donc une méthode d'adaptation locale du maillage, de type hr. Cette approche va dans le sens de la réduction des erreurs de transfert, puisqu'aucun transfert de la solution n'est nécessaire sur les éléments préservés par l'adaptation locale du maillage. De la même façon, cette approche accélère la génération des différents maillages. Par ailleurs, des critères d'arrêt sont utilisés pour déclencher automatiquement les cycles d'adaptation de maillage lorsque cela s'avère nécessaire. Ces critères nécessitent une estimation de l'erreur spatiale à chaque pas de temps, et impliquent la surveillance de plusieurs quantités scalaires dérivant de cette estimation. L'intégration temporelle adaptative repose quant à elle sur un intégrateur temporel hp-adaptatif, basé sur les méthodes BDF. Il dispose de deux mécanismes principaux. En premier lieu, une estimation de l'erreur de troncature locale permet d'ajuster la taille h du pas de temps à l'issue de chaque pas de temps convergé. Dans l'éventualité où la cible renseignée pour cette erreur ne serait pas atteinte, la solution du pas de temps courant peut être rejetée, afin d'intégrer à nouveau le pas de temps avec une valeur réduite de h. Un second mécanisme permet d'ajuster l'ordre p de l'intégrateur temporel (soit le nombre de pas de la méthode BDF) en fonction de considérations fournies par un indicateur de stabilité. L'ordre de l'intégrateur est au demeurant toujours compris en 1 et 3. Ce deuxième mécanisme régule l'utilisation de la méthode d'ordre 3, qui ne possède pas la propriété de stabilité absolue. Les méthodes BDF d'ordre 2 et 3 n'étant pas auto-démarrantes, nous proposons dans ces travaux une approche permettant, postérieurement à chaque adaptation de maillage, de relancer l'intégrateur temporel en conservant les dernières valeurs de h et p utilisées. Pour cela, nous procédons au transfert de la solution associée à plusieurs pas de temps antérieurs. Cette approche assure également une forme de continuité pour les mécanismes adaptifs, qui peuvent ainsi opérer à l'issue du premier pas de temps réalisé après une adaptation du maillage. Le transfert de la solution d'un maillage à l'autre a été étudié pour plusieurs méthodes de transfert : l'interpolation linéaire, l'interpolation consistante, et trois formulations de la projection de Galerkin (ou projection L2). La première formulation de la projection de Galerkin repose sur l'usage d'un maillage auxiliaire (le supermesh), construit par intersection du maillage courant et du maillage résultant de son adaptation. Lorsque la frontière du domaine de calcul est composée exclusivement de segments de droite, cette formulation permet de minimiser l'erreur de transfert en norme L2, en plus de conserver les intégrales sur le domaine des variables dépendantes. La deuxième formulation de la projection de Galerkin appliquedes quadratures de Gauss sur l'élément de référence, directement à partir des éléments du maillage adapté. Bien que globalement moins précise cette deuxième formulation, moyennant l'usage d'un nombre suffisant de points de Gauss, permet néanmoins l'obtention de résultats plus précis à proximité des frontières du domaine (ou des interfaces entre fluides) lorsque celles-ci sont courbes. Partant de ce constat, nous avons également proposé une troisième formulation de la projection, qualifiée d'hybride car issue de la combinaison des deux premières formulations. La deuxième formulation est appliquée sur tous les éléments de bord du maillage adapté, tandis que les éléments restants sont traités via la première formulation. De plus, nous avons implémenté une projection de Galerkin 1D pour post-traiter la solution transférée sur les frontières du domaine de calcul. Nous nous sommes aussi intéressés à une situation particulière : l'usage conjoint d'éléments isoparamétriques d'ordre 2 et d'une formulation ALE des équations de Navier-Stokes. Dans cette situation les éléments du maillage courant, au moment de l'adaptation de maillage, sont délimités par des portions de paraboles. Assurer un transfert de solution suffisamment précis dans ces conditions est une tâche extrêmement ardue. Nous proposons donc dans ces travaux une stratégie visant à contourner en partie cette difficulté. La construction d'un maillage hybride, comportant tant des éléments sub-paramétriques qu'isoparamétriques, est préconisée. Finalement, nous avons développé un transfert local de la solution, afin de tirer avantage de l'adaptation locale du maillage. Cette approche permet une réduction substantielle du coût de calcul lié au transfert de la solution, rendant d'autant plus viable la projection de Galerkin. Faisant suite à une description détaillée des différentes méthodes évoquées, de nombreux résultats de vérification et d'analyse sont également présentés dans cette thèse. Nous avons vérifié en profondeur l'implémentation éléments finis de travail, ainsi que l'implémentation de certains outils développés dans le courant de nos travaux, dont la projection de Galerkin. Plusieurs solutions manufacturées ont été utilisées pour ce faire. Des études de convergence ont été menées pour l'erreur spatiale et temporelle, mais aussi pour l'erreur résultant du transfert d'une solution entre deux maillages. Cette dernière étude a démontré la supériorité de la projection de Galerkin vis-à-vis des méthodes d'interpolation. Nous avons également analysé le comportement des composantes de notre méthode d'adaptation de maillage locale. La formations des poches à remailler, et la relaxation locale du maillage ont notamment été étudiées. Nous illustrons finalement l'utilisation de notre procédure de résolution en considérant quatre écoulements instationnaires différents : un écoulement laminaire autour d'un cylindre circulaire, un jet laminaire impactant sur une paroi chauffée, un écoulement turbulent dans un canal subissant une brusque variation de section (marche descendante), et l'instabilité de Rayleigh-Taylor entre deux fluides séparés par une interface explicite. La variété de ces problèmes, et la richesse de la physique impliquée dans chacun d'eux, nous a permis d'éprouver la robustesse et la pertinence des composantes de notre procédure de résolution. Les résultats montrent que la méthode d'adaptation de maillage locale et les critères d'arrêt employés sont particulièrement efficaces pour assurer la capture très progressive des régions critiques dans l'écoulement, tout en sauvegardant une large partie des éléments du maillage à chaque adaptation. La projection de Galerkin a également permis des transferts de solution de qualité, facilitant ainsi le contrôle du niveau de l'erreur spatiale au cours de la simulation. L'intégrateur temporel adaptif, bien qu'ayant démontré sa force lors de l'usage d'un maillage fixe, a cependant vu son fonctionnement sensiblement perturbé par les adaptations de maillage. Nous avons systématiquement observé une recrudescence de l'erreur de troncature locale lors des reprises du calcul. En conséquence, la taille h du pas de temps s'en est trouvée réduite à l'issue du premier pas de temps, après plusieurs rejets de celui-ci. Les erreurs commises lors du transfert de la solution, mais aussi le fait qu'une solution transférée n'est pas pleinement convergée sur le maillage adapté, en sont la cause. Pour ces mêmes raisons, nous avons constaté que l'indicateur de stabilité verrouillait, en tout temps ou presque, l'accès à la méthode BDF d'ordre 3. Il est à noter que ces limites ne font que brider l'efficacité de notre procédure. La précision des résultats n'est en revanche pas affectée. La réduction automatique de h lors des reprises permet au contraire de limiter considérablement les erreurs de nature itérative. Au demeurant, l'adaptation de la taille h du pas de temps s'est avérée efficace dans l'ensemble. En l'état, notre procédure de résolution peut donc être appliquée pour simuler précisément, efficacement, et automatiquement des écoulements complexes de portée industrielle, mais aussi des écoulements faisant intervenir des mécanismes physiques plus fondamentaux. De plus, les écoulements considérés peuvent être monophasiques, ou polyphasiques à phases séparées.