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.

Abstract

We present in this thesis the development of a general solving procedure for unsteady flows of Newtonian fluids, based on the finite element method. However, as this work is more like a proof of concept, the procedure has been designed in a 2D context only. The overall objective of this procedure is to address the neglected difficulty of simulating such flows both accurately and efficiently. To this end, adaptive methods in space (mesh adaption) and in time (adaptive temporal integration), the relevance of which has been proven in many former works, are used. In doing so, we intend to guarantee the accuracy of the results, and to optimize the computing cost of the simulation. The adaptive procedure thus set up is completely automatic. Indeed, no user intervention is needed during the simulation, as the mesh adaption cycles are triggered automatically. Similar procedures have already been suggested in the literature, but the innovation of ours lies in the choice of its components and in the way they work. We have also emphasized a lot the issues resulting from the transfer of a finite element solution between two meshes. As this operation occurs after each mesh adaption, it is essential to control the errors that arise inevitably. Mesh adaption cycles, performed with a h-adaptive and a r-adaptive methods, allow to locally adjust the mesh size during the simulation, while preserving its overall quality. The h-adaptive method relies on an estimate of the spatial error (committed on each element of the computational domain), derived from a gradient recovery method. Depending on the multiphysics nature of the simulated problem, the spatial error can be translated into one or several error norm(s). Then a transition operator, formulated so to equidistribute the spatial error over the elements of the domain, is used to compute the optimal size of each element of the mesh. Moreover, we used the transition operator to identify the most problematic elements in the mesh, i.e. those whose the optimal size is significantly different from their current size. By removing only these elements and their neighbors from the current mesh, we create pockets (also called cavities) that we can discretize again with elements of a more suitable size. An advancing front meshing tool is used to carry out this last step. The r-adaptive method, which is a combination of several mesh relaxation procedures, is finally applied to maintain the mesh overall quality. Here again, we propose a local formulation of this second method, based on a well-chosen quality metric. The mesh relaxation can thus be carried out only in the regions of the domain that require it the most. Altogether, we are therefore proposing and using a local hr-adaptive mesh adaption method. This approach reduces transfer errors, since a transfer of the solution is only necessary on the elements that do not have an exact match in the previous mesh. This approach also makes the mesh generation more efficient.In addition, stopping criteria are used to automatically trigger mesh adaption cycles when necessary. These criteria imply the monitoring of several scalar quantities derived from an estimate of the spatial error at each time step. Adaptive time integration is carried out with a hp-adaptive time integrator, based on the BDF methods. It has two main mechanisms. First, an estimate of the local truncation error is used to adjust the time step-size h at the end of each converged time step. In the event that this error is too large, the solution of the current time step can be rejected, so it can be computed again with a smaller h value. A second mechanism allows to adjust the time integrator order p (i.e. the number of steps involved in the BDF method). This is done using a stability indicator. Still, the order of the integrator can only vary between 1 and 3. This second mechanism enables a safe use of the method of order 3, kwown to be conditionnaly stable only. Since the BDF methods of order 2 and 3 are not self-starting, we propose in this work an approach to relaunch the time integration efficiently after each mesh adaption. This approch aims at restarting the temporal integration with the last used values of h and p. To do so, we transfer not only the last solution but also those of previous time steps. This approach also makes sure that the adaptive mechanisms can operate immediately after each mesh adaption cycle. The solution transfer from one mesh to another has been studied for several transfer methods : linear interpolation, consistent interpolation, and three formulations of the Galerkin projection (also known as L2 projection). The first formulation of the Galerkin projection makes use of an auxiliary mesh (the supermesh), obtained by intersection of the current mesh and the mesh resulting from its adaption. When the boundary of the computational domain is composed exclusively of line segments, this formulation allows to minimize the transfer error computed in L2 norm, in addition to conserving the integral values of the dependent variables. The second formulation of the Galerkin projection applies Gauss quadrature rules on the reference element, directly from the elements of the adapted mesh. Although globally less accurate, this second formulation, with the use of a sufficient number of Gauss points, nevertheless provides more accurate results close to the domain boundaries (or fluid interfaces) when these are curved. Based on this observation, we have also proposed a third formulation of the projection, called hybrid because it is a combination of the first two : the second one is applied on all the edge elements of the adapted mesh, while the remaining elements are treated with the first one. Besides, we have implemented a 1D Galerkin projection to post-process the transferred solution on the boundaries of the computational domain. We also addressed a particular situation : the use of an ALE formulation of the Navier-Stokes equations together with isoparametric elements of order 2. In this situation elements of the current mesh end up being delimited by portions of parabolas at the time amesh adaption cycle goes off. Ensuring a sufficiently accurate solution transfer under these conditions is an extremely difficult task. We therefore propose in this work a strategy to partially circumvent this difficulty. The construction of a hybrid mesh, including both subparametric and isoparametric elements, is recommended. Finally, we have developed a local transfer of the solution, in order to take advantage of our local mesh adaption method. This approach allows a substantial reduction of the solution transfer computational cost, making the Galerkin projection all the more viable. Many verification and analysis activites where carried out as parts of this research project. Following a detailed description of the different methods already mentioned, this thesis also presents some of the results we thereby obtained. We have verified in depth the finite element implementation, as well as the implementation of some tools developed in the course of our work, the Galerkin projection included. Several manufactured solutions have been desgined to perform this task. Convergence studies were conducted for the spatial and temporal errors, but also for the error arising from the solution transfer between two meshes. The latter study demonstrated the superiority of the Galerkin projection over interpolation methods. We have also analyzed the behavior of the components of our local mesh adaption method. In particular, the formation of the pockets to be remeshed, and the local mesh relaxation were studied. We finally illustrate the use of our solving procedure by considering four different unsteady flows : a laminar flow around a circular cylinder, a laminar impinging jet on a heated wall, a turbulent flow in a channel undergoing a sudden variation of section (backward facing step), and the Rayleigh-Taylor instability between two fluids separated by an explicit interface. The variety of these problems, and the complexity of the physics involved in each of them, allowed us to assess our procedure robustness and relevance. The results show that the local mesh adaption method and the stopping criteria used are particularly effective when it comes to gradually capture the critical regions in the flow, while saving a large part of the mesh elements at each adaption. The Galerkin projection also resulted in accurate solution transfers, making it easier to control the global spatial error during the simulation. The adaptive temporal integrator has however seen its functioning appreciably disturbed by the mesh adaptions.We have systematically observed a recrudescence of the local truncation error each time the computation was resumed. As a result, the time step-size h was automalically reduced following the completion of the first time step, which thus has to be computed again several times. The errors made during the transfer of the solution, but also the fact that a transferred solution is not fully converged on the adapted mesh, are the cause. For these same reasons, we noted that the stability indicator locked, at all times or almost, the access to the BDF method of order 3. It should be noted that these limitations only prevent our solvingprocedure from being even more efficient. As a matter of fact, the decrease of the time step-size h acts as a shield against iterative errors, and thus increases the results accuracy. Moreover, the adaption of the step-size h proved to be quite effective on the whole. As it stands, our solving procedure can therefore be applied to accurately, efficiently and automatically compute complex flows of industrial scope, but also flows involving more fundamental physical mechanisms. In addition, the flows considered can either be monophasic, or multiphasic with segregated phases.