Soutenance mémoire
Télécharger le fichier pdf d’un mémoire de fin d’études
Méthodes de forçage direct pour la simulation des écoule-ments turbulents
La méthode de forçage direct a été introduite par Mohd-Yusof [100–102] bien que ce dernier ne lui ait pas donné cette dénomination. En effet, Fadlun et al. [39] ont été les pre-miers à employer la dénomination direct forcing qui s’est depuis imposée dans la littérature [61, 147, 161]. Les auteurs la qualifient ainsi car, à la différence des approches continues pré-sentées auparavant, la valeur cible est imposée au niveau de la frontière immergée sans passer par un processus de rétro-action. La condition limite est imposée à chaque sous-itération et est ainsi indépendante des fréquences de l’écoulement. Par ailleurs, cette méthode ne souffre pas non plus des limitations du pas de temps mentionnées pour les approches continues ce qui constitue un avantage dans l’optique de son application dans le cadre de simulations instationnaires. Finalement, l’approche par forçage direct est souvent considérée comme la plus à même d’être employée pour des écoulements à haut nombre de Reynolds [90].
En règle général, le terme de forçage direct est déduit à partir d’une discrétisation de la forme incompressible de l’équation du mouvement. Néanmoins, ici, une formulation plus générale du terme source est proposée. Tâchons en premier lieu de considérer une formulation discrétisée générique des équations de Navier-Stokes : [IT]{ΔW} = {R} (1.5).
avec [IT ] l’opérateur d’intégration temporel, W le vecteur des variables conservatives, n l’itération courante, {ΔW } = W n+1 − W n et {R} représente le résidu explicite. L’idée est donc d’imposer W S, le vecteur des variables conservatives cibles, en un point donné du maillage via l’introduction du terme source {F } ({ΔW S} = W S − W n) :
[IT]{ΔW} = {R} + {F} (1.6).
{F} = −{R} + [IT]{ΔWS} (1.7).
En combinant les équations (1.6) et (1.7), on a alors : W n+1 = W S (1.8).
En outre, en prenant [IT ] = Δ1t [I] avec Δt le pas de temps et [I] la matrice identité, W = u avec u le vecteur vitesse et {R} = rhs avec rhs la forme discrétisée des flux explicites, il est possible de se ramener à la formulation plus couramment employée pour le forçage direct [61] : u n+1 − u n = rhs + f (1.9) Δt u S − u n f = − rhs + (1.10).
Au sein même de la méthode de forçage directe, une distinction apparaît selon la nature directe ou indirecte de la technique d’imposition des conditions limites. Dans sa formulation d’origine [100, 101] le terme source proposé apparaît comme un terme à calculer explicite-ment (indirect BC imposition) [71, 153, 161, 172]. Dans la plupart des implémentations de la méthode, le terme source n’est cependant jamais calculée explicitement mais son effet (cf. equation (1.8)) est pris en compte de manière implicite en reconstruisant les variables de vitesse au niveau de certains points au voisinage de la paroi [147] (direct BC imposition). Une telle distinction engendre donc des stratégies de représentation de la paroi pouvant être radicalement différentes comme cela est explicité dans le paragraphe suivant.
Prise en compte de la paroi au sein du maillage de fond : L’objectif d’une méthode de frontière immergée est de reproduire les effets des parois d’un solide sur l’écoulement. Pour ce faire, les termes sources mentionnés précédemment (les méthodes de maillage car-tésien ne sont pas concernées) permettent d’imposer la vitesse des parois du solide sur les points du maillage de fond. Cependant, en pratique, très peu de points de la grille de fond coïncident avec la paroi modélisée (cf. figure 1.2(b)). Le type de forçage le plus simple consiste alors à faire uniquement la distinction entre les points à l’intérieur et à l’extérieur du domaine Ωs sans prendre en compte la position exacte de la frontière . Il découle de cette opération une légère transformation de la géométrie réelle considérée car les contours prennent la forme du maillage environnant. Dans ce cas, on parle de description de la géo-métrie en marche d’escalier (stepwise) [39]. Avec cette approximation, la position exacte de la frontière modélisée n’est connue qu’à une maille près et engendre une précision du calcul inférieure à l’ordre un [39]. À cet égard, Lamballais [78] insiste sur le fait que si aucun soin n’est pris pour prendre en compte la position exacte des frontières, aucune des méthodes de modélisation de paroi par forçage ne peut prétendre à un ordre de précision supérieur
à un, sauf naturellement dans les cas très particuliers où les parois considérées sont situées exactement sur les points de la grille. Il est à noter que ce manque de précision ne concerne 42 efforts aérodynamiques dans le contexte des techniques de frontières immergées en pratique que la région de très proche paroi.
À l’image des méthodes de frontières immergées, les solutions proposées pour recons-truire la position de la paroi sont nombreuses et peuvent être divisées en deux catégories principales [61, 147]. La première regroupe les méthodes qui affectent une pondération au terme source dans le voisinage de la frontière immergée. En général, plusieurs nœuds de maillage sont concernés au voisinage de la paroi dont la représentation est alors qualifiée de diffuse (diffuse interface methods). Dans la seconde catégorie, l’opération n’est pas effec-tuée sur le terme source mais directement sur la variable que l’on cherche à reconstruire. Ce type de méthode fait appel à des interpolations permettant une représentation précise de la paroi (sharp interface methods).
Application à la simulation numérique des écoule-ments turbulents
Cette section vise à présenter les possibilités offertes par les frontières immergées pour traiter les écoulements turbulents à hauts nombres de Reynolds. S’inspirant de Iaccarino et Verzicco [61], la première partie permet d’illustrer le développement historique des appli-cations faisant intervenir des IBC. Dans un seconds temps, les problématiques spécifiques qu’elles engendrent sont plus particulièrement abordées ainsi que les solutions proposées dans la littérature.
Évolution historique des simulations d’objets rigides immer-gés
Comme cela a déjà été mentionné dans l’introduction de la section précédente, la formu-lation originelle des frontières immergées a été introduite par Peskin [118, 119] au début des années 70 afin de décrire l’écoulement sanguin dans le cœur. Il s’agit d’une méthode conçue pour tenir compte des déformations élastiques d’une paroi dans un écoulement incompres-sible à de faibles nombres de Reynolds. En principe, cette approche peut être appliquée à une paroi solide en réduisant le paramètre modélisant son caractère déformable. En pratique, cependant, une telle opération mène à des instabilités numériques [152]. En conséquence, elle a depuis été principalement employée pour des applications tournées vers la Bioméca-nique [6, 41, 121]. Il faut attendre la décennie suivante pour voir apparaître les premières applications de l’approche à des parois solides et indéformables [61]. Briscolini et Santan-gelo [11] emploient une approche de frontières immergées (Mask Method) pour calculer des écoulements bidimensionnels instationnaires autour d’obstacles cylindriques et carrés pour des nombres de Reynolds allant jusqu’à 1000.
Le début des années 90 voit l’introduction d’une nouvelle méthode de frontières immer-gées spécifiquement dédiée à la prise en charge de parois solides et indéformables. Il s’agit de la méthode des frontières virtuelles de Goldstein et al. [51], brièvement décrite dans la partie précédente. Ces derniers ont alors simulé l’écoulement bidimensionnel généré par le démarrage d’un cylindre et l’écoulement turbulent passant au sein d’un canal aux parois lisses d’une part et striées d’autre part. Les simulations sont réalisées dans le cadre d’approches spectrales et le terme source est appliqué au moyen d’une fonction de lissage. Cette fonction est introduite pour réduire l’intensité des oscillations numériques parasites générées à la paroi. Saiki et Biringen [136] ont par la suite employé la même fonction de forçage pour simuler l’écoulement autour d’un cylindre avec et sans mouvement de rotation. L’approche numérique qu’ils utilisent est une méthode aux différences finies centrées d’ordre quatre. Un tel schéma permet d’atténuer l’apparition des oscillations numériques proche de la paroi. Les auteurs emploient par ailleurs une interpolation bilinéaire afin de déterminer la vitesse pour un point appartenant à la paroi immergée à partir des quatre points du maillage de fond environnants (en 2D). Ils qualifient ce schéma d’interpolation de précis à l’ordre un et similaire à la fonction delta de Peskin [120]. Toutefois, en raison des inconvénients énoncés dans la première partie (présence de deux constantes dépendantes du problème considéré, limitation du pas de temps pour des écoulements instationnaires et réduction de la précision en raison de la fonction de lissage) le nombre de problèmes traités avec cette approche est resté relativement restreint [61].
A la fin des années 90, Mohd-Yusof [101] élabore une nouvelle formulation du terme source qui permet de lever les difficultés posées par les formulations antérieures. La fonction de forçage introduite, qualifiée par la suite de forçage direct [39], n’affecte pas la stabilité du schéma d’intégration temporelle employée et ne nécessite pas l’emploi d’une fonction de lissage. Par ailleurs, cette formulation n’introduit aucune constante dépendante de l’écoule-ment. Dans son étude, Mohd-Yusof associe l’introduction de cette fonction de forçage avec l’emploi de B-splines pour simuler l’écoulement laminaire à l’intérieur d’un canal au paroi striée. Les avantages tirés de cette formulation sont remarquables vis-à-vis des précédentes formulations [61]. Les équations de Navier-Stokes incompressibles sont résolues au moyen d’une méthode spectrale et ont été appliquées sur des simulations d’écoulements autour de cylindres et de sphères pour des nombres de Reynolds modérés [61]. L’idée a par la suite été mise en œuvre au début des années 2000 par Fadlun et al. [39] et Verzicco et al. [165] au sein d’un code de différences finies et appliquée à la simulation des grandes échelles (LES) de géométries simples et complexes pour une large gamme de nombres de Reynolds. Ces études ont permis de mettre en avant les attraits de la méthode de forçage direct, notam-ment lorsqu’elle est employée dans le cadre de la LES.
Au cours des années 2000, le nombre d’études impliquant la technique des frontières immergées connaît un essor considérable. Ces méthodes s’imposent comme un moyen po-lyvalent de traiter des géométries très complexes [98], comme en témoigne par exemple la simulation du vol d’une libellule (cf. figure 1.11). Aujourd’hui, les études impliquant une technique de frontière immergée couvrent un large spectre d’applications [70] : écoulements compressibles [45, 91], écoulements de particules [161, 191], microfluidique [3], interaction avec des corps solides ([48, 194]), écoulements multiphasiques [44], prise en compte des trans-ferts thermiques [60, 192], écoulements environnementaux [146], fluides biologiques [40], etc. Néanmoins, en dépit du nombre très élevé d’études (au regard des décennies ultérieures) ayant été publiées dans le courant des années 2000, seul un nombre relativement restreint d’entre elles traitent de la thématique des écoulements turbulents pour des hauts nombres de Reynolds [70, 196].
La simulation des écoulements pour des hauts nombres de Reynolds nécessite l’intro-duction d’un dispositif de modélisation de la turbulence. La mise en place de frontières immergées dans un tel contexte implique donc un traitement particulier des variables tur-bulentes afin de ne pas détériorer la précision et la stabilité du calcul [196]. Dès 2002, Kalitzin et Iaccarino [67] emploient la méthode de forçage direct implicite de Fadlun et al. avec les modèles de turbulences de Spalart-Allmaras [149] et de Wilcox [179]. Ils réalisent la simula-tion d’un écoulement autour d’une aube de turbine pour un nombre de Reynolds Re = 148 000. Cependant, cette technique n’est valable que si le maillage est suffisamment raffiné au niveau de la paroi, ce qui n’est pas forcément le cas lorsque des frontières immergées sont introduites dans une grille de fond. Des méthodes permettant de pallier ce problème ont par la suite été proposées et sont exposées dans la section suivante. La simulation des grandes échelles est une autre option pour simuler la turbulence. Tou-tefois, au sein de la thématique des frontières immergées, la simulation des grandes échelles semble occuper une place plus réduite avec un nombre de publications estimé à 10% 1 de l’ensemble des études relatives aux IBC. Et même si l’engouement pour les techniques de simulations avancées de type LES semble s’être accru au cours de la dernière décennie, le nombre d’études traitant des écoulements turbulents à haut nombre de Reynolds est resté limité [70]. Ainsi, Tessicini et al. [155] ont simulé l’écoulement autour du bord de fuite d’un hydroglisseur pour un nombre de Reynolds basé sur la vitesse à l’infini amont et la corde c du profil Rec = 2, 15 × 106 . Ce calcul a mis en jeu des frontières immergées avec une technique LES impliquant un modèle de loi de paroi assez simple [13, 170, 171]. Eisenbach et Friedrich [38] ont par ailleurs simulé l’écoulement autour d’une aile pour un nombre de Reynolds Rec = 105 et un angle d’incidence élevé α = 18◦ (cf. figure 1.12). En dépit de ces applications des difficultés subsistent pour ce type de simulations et sont également abordées dans la section suivante.
Finalement, les études traitant des frontières immergées dans le cadre de simulations avancées hybrides RANS/LES sont très rares [46]. En particulier, aucune ne traite de l’ap-plication des frontières immergées associée à la Zonal Detached Eddy Simulation (ZDES) [29, 30].
Simulation des écoulements turbulents à haut nombre de Reynolds
Le principal attrait d’une approche de frontières immergées réside dans sa capacité à prendre en compte une géométrie complexe sans introduire les contraintes propres aux maillages adaptés à la forme de l’objet. L’effet de la paroi sur l’écoulement est reproduit via l’action locale d’un terme source ou la reconstruction directe de la solution à proxi-mité de la paroi. L’indépendance du maillage vis-à-vis de la forme de l’objet considéré offre notamment la possibilité de tirer profit de la robustesse et de la précision des solveurs struc-turés. L’avantage consiste donc à pouvoir réduire le temps dévolu à l’étape de génération du maillage pouvant se révéler à la fois coûteuse et fastidieuse dans le cadre de maillages adaptés à la forme de l’objet. Toutefois, la conséquence directe de ce découplage entre la paroi et le maillage est une difficulté à assurer une résolution satisfaisante du maillage dans les régions à forts gradients de vitesse à proximité d’une paroi solide. En effet, contrairement aux maillages adaptés à la forme de l’objet qui permettent facilement d’accroître le nombre de points dans une direction donnée uniquement, l’augmentation de la résolution engendrée par la proximité d’une paroi s’opère souvent dans deux ou trois directions lorsqu’il s’agit d’une méthode de frontières immergées. Dans le cadre des écoulements à hauts nombres de Reynolds, le nombre de points générés pour pouvoir résoudre correctement une couche limite turbulente peut très vite devenir prohibitif en terme de coût de calcul.
Parmi les solutions proposées figurent celles qui ont trait à la nature du maillage de fond. La première option consiste à employer une stratégie de raffinement du maillage par adaptation locale (Adaptive Mesh Refinement (AMR) strategy) à proximité de la paroi. Cette solution existe plus généralement pour les maillages structurés [66]. Elle a été introduite dans le contexte des IBC par Kalitzin et Iaccarino [68] et reprise par la suite [14, 52, 164]. Elle a notamment été employée par Kang et al. [70] dans le cadre d’une si-mulation d’un écoulement autour d’un profil d’aile pour un nombre de Reynolds basé sur la vitesse à l’infini amont et la corde Rec = 150 000 et un angle d’attaque égal à 8◦ . Dans ce cas, la technique de raffinement local permet de réduire de 70% le nombre de points nécessaires par rapport à un maillage purement cartésien (cf. figure 1.13). Cet avantage est néanmoins obtenu au prix d’une procédure de gestion du maillage non structuré (octree) pouvant tout de même se révéler lourde et coûteuse à mettre en œuvre [128].
Prévision des efforts aérodynamiques en présence de parois immergées
Dans un cadre industriel, la simulation numérique des écoulements autour d’un objet a souvent pour finalité d’évaluer les efforts aérodynamiques qui s’exercent sur ce dernier. Dans le contexte applicatif qui nous intéresse, le chargement dynamique joue un rôle prépondérant. Cette section vise donc à s’interroger sur le comportement de la pression au voisinage d’une paroi immergée ainsi que sur la manière de reconstituer ces efforts en pratique.
Comportement et traitement de la pression à proximité de la paroi
Le comportement de la pression au voisinage d’une paroi (et plus généralement d’une interface) immergée est problématique à plus d’un titre. Historiquement, c’est dans le cadre de l’approche continue appliquée à des parois mobiles et déformables que les premières dif-ficultés sont apparues. En effet, Peskin [119] constate que l’application directe de son terme de forçage engendre l’apparition d’ondes de pression parasites prenant naissance au niveau de la paroi immergée. Ces ondes sont provoquées par la discontinuité du champ qui sur-vient au moment du changement d’état (fluide/solide) d’un nœud au passage de l’interface. Dans le cas de la méthode des frontières immergées de Peskin, et plus globalement pour les approches continues, l’introduction d’une fonction de lissage (fonction delta discrète) permet d’atténuer l’apparition de ces ondes numériques. Un tel problème n’est cependant pas propres aux approches continues. En effet, les approches discrètes y sont également confrontées [90]. Concernant les méthodes discrètes qui représentent précisément la paroi, Borazjani et Sotiropoulos [9] ont montré qu’il est possible de réduire considérablement l’am-plitude des oscillations parasites en raffinant le maillage proche de la paroi. En pratique, un tel raffinement n’est pas toujours réalisable et un traitement numérique propre à la méthode de frontière immergée et au schéma de discrétisation est alors nécessaire. Yang et Balaras [186] ont par exemple proposé une approche consistant à étendre les champs de pression et de vitesse dans la phase solide afin d’améliorer la précision des dérivées au voisinage de la paroi. Une version simplifiée de cette approche a été développée dans [187], où seul le gradient de pression est reconstruit pour les nœuds solides proches de la paroi immergée. Par ailleurs, une approche alternative a été développée par Seo et Mittal [140] qui com-binent un forçage direct implicite avec une approche de type cut-cell. Concernant l’origine numérique des ondes parasites, Lee et al. [81] mettent en avant deux causes. La première est une discontinuité de pression provoquée par le passage d’un nœud de l’état solide à l’état fluide. La seconde est engendrée par une discontinuité temporelle sur la vitesse au niveau d’un nœud passant de l’état fluide à l’état solide. Ces discontinuités engendrent alors une perte de la conservation de la masse proche de la paroi [83, 140].
La problématique de la conservation de la masse n’est pas propre aux parois mobiles et/ou déformables mais concerne également les parois fixes et indéformables [70, 147]. Dans ce cas, le non respect de l’équation de continuité se traduit par une dégradation de la précision pour la pression au voisinage de la paroi. Plusieurs solutions à ce problème ont été proposées dans la littérature. Avec une approche aux volumes finis, Kim et al. [72] ont introduit un terme source (ou puits selon le cas) au sein de l’équation de continuité afin de rétablir la conservation de la masse. Une méthode alternative a été proposée par Li et Wang [87] pour un schéma de différences finies. Ces derniers ont trouvé que l’imposition d’un gradient de pression nul dans la direction normale à la paroi permettait d’avoir une distribution correcte de la pression et de respecter la conservation de la masse dans le cadre d’écoulements incompressibles autour de corps non profilés. Par ailleurs, dans le cadre d’une méthode de forçage directe pour les simulations LES/DNS, Kang et al. [70] ont proposé une méthode appelée immersed boundary-approximated domain method. Dans cette méthode (cf. figure 1.19(c)), la conservation de la masse est satisfaite pour l’ensemble des cellules fluides non intersectées par la paroi immergée (zone grisée). Pour ce faire, la vitesse est reconstruite aux bornes de cette zone (frontière noire en marche d’escalier) par une méthode d’interpolation linéaire et un terme de correction est ajouté de sorte à satisfaire la contrainte de conservation de la masse dans la région bornée d’un côté par la région fluide et de l’autre par le corps immergé (représenté par IB). Cette méthode a notamment été employée pour simuler l’écoulement autour d’un profil d’aile par une approche LES. Le nombre de Reynolds basé sur la corde et la vitesse à l’infini amont est Re = 150 000. Les résultats, qui incluent en particulier un spectre de puissance de la pression à la paroi, sont en accord avec des données expérimentales et une autre simulation réalisée à partir d’un maillage adapté à la forme de l’objet. Ces résultats suggèrent que la correction introduite permet de prévoir correctement la dynamique de proche paroi des écoulements turbulents. Finalement, on notera qu’au sein des méthodes de frontières immergées au sens large, les méthodes de type cut-cell (cf. figure 1.19(b)) ne sont pas concernées par la problématique de conservation de la masse. Il s’agit cependant d’un type de méthode très difficile à mettre en œuvre sur des cas tridimensionnels et qui impose, pour être intégré à un code de calcul structuré, de lourdes modifications de ce dernier.
Reconstruction des efforts aérodynamiques
Le calcul des efforts aérodynamiques est une procédure de post-traitement courante en CFD. Dans le cadre de maillages adaptés à la forme de l’objet d’étude, elle est relativement immédiate à mettre en œuvre. En effet, dans un tel contexte, les données à intégrer pour reconstituer la résultante des efforts qui s’exercent sur l’objet peuvent être facilement déter-minées à partir du champ aérodynamique et de la géométrie associés aux points du maillage représentatif du domaine d’étude. En revanche, avec l’emploi d’une méthode de frontières immergées cette procédure est moins triviale en raison de l’absence de lien a priori entre les variables de l’écoulement issues du maillage de fond et les données géométriques relatives à la surface de l’objet immergé. Bien que la thématique des frontières immergées soit large-ment traitée dans la littérature, une telle problématique est rarement abordée [4]. Parmi les techniques proposées, on distingue les approches qualifiées ici de directes [169, 177, 178] et indirectes [4, 21, 75].
L’approche directe consiste à interpoler les variables d’intérêt du champ aérodyna-mique sur les surfaces élémentaires représentatives de l’objet immergé. Dans ce cas, deux problématiques se posent, à savoir la détermination des surfaces élémentaires considérées et l’évaluation des variables d’intérêt au niveau de ces surfaces. Dans le cadre de problèmes d’interaction fluide-structure, Wang et al. [169] ont proposé deux algorithmes de reconstruc-tion du chargement dynamique qui s’exerce sur une surface immergée. Dans ce cas, l’objectif des auteurs est d’avoir une communication des conditions limites entre le code de mécanique des fluides et celui de mécanique des solides employés. Bien qu’il ne s’agisse pas de l’objectif visé dans cette thèse, les idées développées semblent tout de même pouvoir s’appliquer dans le contexte de ces travaux. Sans entrer dans les détails, le premier algorithme présenté a pour tâche de redéfinir les éléments du maillage d’origine afin que leur taille soit commensu-rable à la taille des cellules du maillage de fond (cf. figure 1.20 à gauche). Les variables sont alors déterminées sur ces éléments en tant que solutions d’un problème de Riemann à une dimension. L’autre proposition (cf. figure 1.20 à droite) consiste à approximer la paroi par sa représentation en marche d’escalier. Les variables de l’écoulement sont, dans ce cas, di-rectement issues du maillage de fond. Parmi ces deux solutions, la première semble la plus à même de reconstituer avec précision les charges latérales. Néanmoins, c’est également celle qui apparaît la plus délicate à mettre en œuvre, particulièrement pour des cas tridimen-sionnels. En effet, la mise en place d’un algorithme de subdivision du maillage surfacique d’origine permettant d’accorder la taille des triangles représentatifs de la paroi avec la taille des cellules du maillage de fond, bien que possible, semble demander un effort important. Il s’agit alors de savoir si l’approximation engendrée permet tout de même de reconstituer avec précision les charges latérales dans le cadre de l’écoulement qui nous intéresse dans ces travaux.
Présentation générale du code RAYTRACER3D
À des fins d’appropriation de la méthode et pour maximiser la compatibilité des entrées et sorties avec le code de calcul FLU3M, des algorithmes de tracé de rayon et de calcul des distances à la paroi ont été développés dans un code nommé RAYTRACER3D. Le développement de cet outil de prétraitement au cours de la thèse a permis un libre choix de la nature du terme source à imposer et surtout de son lieu d’application. Le principe de fonctionnement du code RAYTRACER3D (illustré par la figure 2.5) se base sur la lecture d’un fichier texte. Ce fichier texte, qui porte le nom générique ibdon, est constitué d’un ensemble de lignes constituant la carte des commandes permettant in fine de réaliser les opérations souhaitées telles que la différenciation fluide/solide du maillage de fond ou le calcul des distances à la paroi.
Pour une utilisation classique du code RAYTRACER3D, la carte présente une succession de commandes qu’il est possible de regrouper de la manière suivante (cf. figure 2.5) :
1. Déclaration générale des bases de données.
L’utilisateur indique au programme combien de domaines, de surfaces à immerger et éventuellement de conditions aux limites classiques vont intervenir lors de l’exécution du code. Cela permet au programme de créer les bases de données employées lors des étapes suivantes pour gérer l’ensemble de la configuration traitée.
Gestion des domaines
FLU3M est un code aux volumes finis structuré et multi-domaines. Le maillage de fond est ainsi représenté par plusieurs domaines reliés entre eux par des conditions de raccord informatique. L’étape qui suit la déclaration des bases de données est la déclaration des dimensions du maillage et le chargement de ce dernier. Un maillage structuré est constitué d’hexaèdres que l’on désignera également par le terme cellule. Chaque cellule est composée de son centre (center) et de ses huit sommets (vertex). L’information des points du maillage est donnée en vertex, alors que l’on calcule les in-formations nécessaires en center. Aussi faut-il procéder à une interpolation du maillage après chargement. Une fois le maillage chargé, on initialise également les distances à la paroi par celles de la configuration de référence (i.e. sans objet immergé)
Gestion des objets à immerger
Cette étape consiste à créer les bases de données propres à un objet à immerger et à charger son maillage surfacique. On rappelle que le maillage surfacique est une collection de triangles non connectés entre eux. Le format STL permet de connaître les coordonnées cartésiennes des trois sommets ainsi que la normale au plan de chaque triangle qu’il est possible de recalculer selon l’orientation souhaitée.
Différenciation fluide/solide
La différenciation fluide/solide est la procédure permettant d’identifier les centres des cellules du maillage de fond qui se trouvent à l’intérieur du corps à immerger. Cette information est ensuite utilisée par le code FLU3M afin d’appliquer le terme source permettant de reproduire l’effet des parois du corps sur l’écoulement (cf. section 2.2.5). Le résultat de la procédure d’identification est alors stocké dans une variable tagibc qui vaut 1 si le centre de la cellule est à l’intérieur du corps et 0 s’il se situe à l’extérieur. Cette étape fondamentale, commune à toutes les méthodes de modélisation de paroi par forçage, fait l’objet d’une description détaillée dans la section suivante 2.1.3.
Calcul des distances à la paroi
Outre la valeur de la variable tagibc, l’application de la modélisation de paroi par forçage nécessite la mise à jour des distances à la paroi du fait de la présence des objets à immerger. Ces distances sont nécessaires pour l’utilisation du modèle de turbulence de Spalart-Allmaras et de la ZDES (cf. section 2.2.2). Le choix d’une méthode de calcul des distances à la paroi dans le contexte de la mise en œuvre des frontières immergées et de la ZDES n’est pas trivial. Il fait donc l’objet d’une description détaillée dans la section 2.1.4.
Extraction des données
Une fois les variables cibles calculées, il est nécessaire de les écrire dans des fichiers afin qu’elles puissent être lues et utilisées par le solveur FLU3M. Cette étape marque la fin de l’opération de prétraitement.
Description et implémentation de la technique du tracé de rayon
L’utilisation d’une méthode de frontières immergées implique que les équations de Navier-Stokes sont résolues sur un maillage dont les points ne coïncident en général pas avec les parois de l’objet dont on cherche à reproduire l’effet. L’étape préalable, commune à la ma-jorité de ces méthodes, consiste alors à identifier les points du maillage de fond qui sont à l’intérieur du solide immergé (points solides), ceux qui sont à l’extérieur (points fluides) et éventuellement les points frontières (points interfaces), comme illustré par la figure 2.6.
Considérations sur l’optimisation du tracé de rayon
L’objectif principal de cette thèse n’était pas d’implémenter un code de tracé de rayon optimisé. Toutefois, on se propose ici d’exposer quelques pistes d’amélioration. La technique la plus couramment employée est certainement celle qui consiste à définir la plus petite boîte cartésienne (bounding box) contenant l’objet à immerger [9, 26]. L’objectif est de réduire le nombre de points sur lesquels l’opération de tracé de rayon est effectivement réalisée. Pour ce faire, les points sont au préalable identifiés par rapport à leur position relativement à cette boîte. Si le point est à l’extérieur, alors on sait qu’il est fluide. Sinon, on réalise la procédure du tracé de rayon. Une telle technique a été implémentée. L’autre technique répandue (Control cell method ou Grid cell method) [9, 26] consiste à subdiviser le domaine d’investigation en boîtes cartésiennes et à commencer par déterminer si ces dernières contiennent des triangles de la géométrie surfacique. La procédure du tracé de rayon est alors appliquée sur chaque point de ces boîtes, les autres étant étiquetés suivant le statut de la boîte qui les contient. Bien qu’une telle procédure semble être à même de réduire drastiquement le temps de calcul, cette dernière n’a pas été implémentée au regard des développements engendrés par sa mise en place.
Description et implémentation du calcul des distances à la paroi
Certains modèles de turbulence requièrent la connaissance des distances à la paroi pour pouvoir être appliqués. C’est le cas du modèle de Spalart-Allmaras [149] utilisé dans cette étude. Ainsi, outre l’étape d’identification des cellules fluides et solides des cellules du maillage de fond, il convient également de leur attribuer les distances à la paroi qui tiennent compte de la présence de l’objet immergé. Rappelons que ce dernier est décrit par un maillage surfacique composé d’une série de triangles sans lien les uns avec les autres. Pour un point donné, il s’agit donc de déterminer l’ensemble des distances entre ce point et les triangles de la surface et de ne conserver que la valeur minimale. Dans cette partie sont donc abordées la méthode de calcul d’une distance dans l’espace entre un point et un triangle, puis le calcul des distances à la paroi pour l’ensemble de la configuration.
Distance minimale entre un point et un triangle
Le calcul de la distance minimale entre un point et un triangle n’est pas immédiat car trois configurations sont possibles suivant la position dans l’espace de ce point. En effet, comme l’illustre la figure 2.11, la distance la plus courte du point P au triangle (Vi)i=1,3 peut aussi bien être la distance au plan du triangle que la distance à l’une de ses arêtes ou que la distance à l’un de ses trois sommets.
Il s’avère qu’un tel problème est fréquemment traité dans le domaine de la représenta-tion graphique par ordinateur [65] où les distances à la paroi sont utilisées afin de faciliter les représentations bidimensionnelles (e.g. un écran d’ordinateur) d’objets tridimensionnels. Dès lors, la littérature propose de nombreux moyens de traiter le problème. Par exemple, Jones [64] présente deux méthodes (méthode 3D et méthode 2D) pour réaliser ce calcul. Bien que la méthode 2D apparaisse être plus rapide d’un facteur quatre que la méthode 3D, c’est cette dernière qui a tout de même été implémentée au cours de la thèse car sa philosophie se rapproche plus de la technique du tracé de rayon (et est ainsi plus facile à mettre en œuvre). Soit P un point quelconque de l’espace et (Vi)i=1,3 un triangle, la technique (illustrée dans la figure 2.12 et formalisée par les algorithmes 4 et 5) permettant de calculer la distance minimale entre P et (Vi)i=1,3 consiste en les étapes suivantes :
1. Déterminer le point P 0 par la projection du point P sur le plan du triangle (Vi)i=1,3.
2. Calculer f1, f2 et f3 avec fi = (OVi ∧ ViP 0) · n, avec O un point strictement à l’intérieur du triangle (par exemple le centre de gravité) et n le vecteur normal au plan du triangle.
3. Déterminer le secteur angulaire VidOVj contenant P 0 (parmi V1dOV2, V2dOV3 et V3dOV1) avec les signes de f1, f2 et f3.
– Si f1 ≥ 0 et f2 ≤ 0 alors P 0 se trouve dans le secteur angulaire V1dOV2.
– Si f2 ≥ 0 et f3 ≤ 0 alors P 0 se trouve dans le secteur angulaire V2dOV3.
– Si f3 ≥ 0 et f1 ≤ 0 alors P 0 se trouve dans le secteur angulaire V3dOV1.
4. Déterminer si P 0 est à l’intérieur ou à l’extérieur du triangle VidOVj.
– Si gij = P 0Vi ∧ P 0Vj · n < 0, P 0 est à l’extérieur du triangle P.
– est à l’intérieur et la distance minimale entre le point et le triangle. Sinon ilest ||P P 0|| (distance point-plan).
5. Déterminer P 00 projection de P 0 sur la droite (ViVj) et calculer t tel que ViP 00 = tViVj.
– Si 0 ≤ t ≤ 1, P 00 est sur le segment ViVj et la distance minimale entre le point P et le triangle est ||P P 00||.
– Si t < 0, P est plus proche de Vi et la distance minimale entre le point P et le triangle est ||P Vi|| (distance point-sommet).
– Si t > 1, P est plus proche de Vj et la distance minimale entre le point P et le triangle est ||P Vj|| (distance point-sommet)(distancepoint-arête).
Modélisation de la turbulence : Zonal Detached Eddy Simulation (ZDES)
Du fait de leur non linéarité, la résolution des équations de Navier-Stokes n’est possible analytiquement que dans certains cas particuliers et sous certaines hypothèses fortes. Afin d’obtenir une représentation de l’écoulement considéré, il est possible d’avoir recours à la simulation numérique. Les écoulements turbulents sont caractérisés par une large gamme d’échelles spatiales et temporelles dont la simulation numérique directe (DNS pour Direct Numerical Simulation) reste, aujourd’hui encore, inenvisageable pour la plupart des cas d’application en aéronautique en termes de coût de calcul et ce en dépit de l’avènement des supercalculateurs et de l’augmentation de la puissance des moyens de calcul. En ef-fet, d’après Sagaut et al. [135] le temps de calcul de la résolution directe d’un problème tridimensionnel est estimé à tcalcul ∝ Re 3L (avec ReL le nombre de Reynolds basé sur L une longueur caractéristique du problème traité). Pour des écoulements à haut nombre de Reynolds, une telle estimation implique de devoir recourir à la modélisation de la turbulence.
La modélisation de la turbulence repose sur une séparation des échelles de la turbulence qui résulte d’un compromis entre les échelles résolues les plus susceptibles d’influencer la dynamique de l’écoulement et les échelles non résolues dont l’influence est prise en compte par un modèle mathématique. On peut ainsi distinguer cinq grandes classes de méthodes suivant le degré de modélisation choisi (cf. figure 2.17). Suivant le type d’étude effectuée, la méthode employée résulte du compromis entre coût et degré de modélisation. Déjà men-tionnée plus haut, l’approche de type DNS est dépourvue de modèle mathématique pour la turbulence. C’est ainsi l’approche la plus coûteuse. À l’opposé se trouvent les approches qualifiées de RANS (Reynolds Averaged Navier Stokes) qui résolvent les équations de Navier-Stokes moyennées en temps. Dans ce cas, le champ turbulent est entièrement modélisé. Le succès de ce type d’approche réside dans leur capacité à prévoir le champ aérodynamique moyen avec une précision suffisante pour un grand nombre d’applications et ce pour un coût de simulation en accord avec les contraintes industrielles de temps de restitution. De nos jours, l’approche RANS, ainsi que son extension instationnaire (URANS), reste massivement utilisée pour les phases de conception ( e.g. études d’optimisation, couplages multidiscipli-naires, etc). Il existe cependant des cas d’applications, avec la présence de phénomènes de recirculations complexes et fortement instationnaires, où l’approche montre ses limites.
|
Table des matières
1 Simulations numériques des écoulements turbulents et évaluation des efforts aérodynamiques dans le contexte des techniques de frontières immergées
1.1 Introduction de la thématique des frontières immergées
1.1.1 Philosophie, contexte applicatif et typologie des méthodes
1.1.2 Méthodes de forçage direct pour la simulation des écoulements turbulents
1.2 Application à la simulation numérique des écoulements turbulents
1.2.1 Évolution historique des simulations d’objets rigides immergés
1.2.2 Simulation des écoulements turbulents à haut nombre de Reynolds
1.3 Prévision des efforts aérodynamiques en présence de parois immergées
1.3.1 Comportement et traitement de la pression à proximité de la paroi
1.3.2 Reconstruction des efforts aérodynamiques
1.4 Synthèse
2 Définition d’une stratégie numérique zonale couplée à une méthode de frontières immergées
2.1 Description et implémentation de la procédure d’immersion
2.1.1 Description surfacique de l’objet à immerger
2.1.2 Présentation générale du code RAYTRACER3D
2.1.3 Description et implémentation de la technique du tracé de rayon
2.1.3.1 Intersection rayon/triangle
2.1.3.2 Problématiques et modifications de l’algorithme
2.1.3.3 Considérations sur l’optimisation du tracé de rayon
2.1.4 Description et implémentation du calcul des distances à la paroi
2.1.4.1 Distance minimale entre un point et un triangle
2.1.4.2 Distances à la paroi de l’ensemble de la configuration
2.1.5 Exemples de mise en oeuvre
2.2 Mise en oeuvre de la modélisation de paroi par forçage direct au sein de la simulation numérique des écoulements turbulents
2.2.1 Lois de conservation des écoulements turbulents compressibles
2.2.1.1 Équations du mouvement
2.2.1.2 Fermeture du système
2.2.2 Modélisation de la turbulence : Zonal Detached Eddy Simulation (ZDES)
2.2.2.1 Motivations
2.2.2.2 Modèle de turbulence de Spalart-Allmaras
2.2.2.3 Detached Eddy Simulation
2.2.2.4 Zonal Detached Eddy Simulation
2.2.3 Description générale du code de calcul FLU3M
2.2.3.1 Intégration en temps
2.2.3.2 Flux de convection et de diffusion
2.2.3.3 Traitement de la phase implicite
2.2.4 Traitement des conditions de paroi classiques
2.2.4.1 Initialisation des cellules fictives pour le calcul du flux convectif
2.2.4.2 Calcul du gradient dans les cellules fictives pour le calcul des flux visqueux
2.2.5 Implémentation de la technique de forçage direct au sein de FLU3M .
2.2.5.1 Définition du terme de forçage direct
2.2.5.2 Imposition des conditions limites 2.3 Synthèse
3 Étude des aptitudes de la méthode à reproduire les effets de la paroi et application à des configurations de lanceur
3.1 Études numériques préliminaires de la formulation IBC
3.1.1 Écoulement de couche limite turbulente sur une plaque plane immergée103
3.1.1.1 Présentation du cas test numérique
3.1.1.2 Mise en oeuvre de la méthodologie
3.1.1.3 Étude de l’effet du terme source dans la région de proche paroi106
3.1.2 Écoulement de sillage derrière une plaque plane immergée
3.1.2.1 Présentation et mise en oeuvre du cas test numérique
3.1.2.2 Étude des capacités du terme source à restituer un effet de blocage
3.1.3 Synthèse des deux études numériques préliminaires
3.2 Application et évaluation de la stratégie zonale sur une configuration simplifiée de lanceur : immersion d’une jupe dentelée cylindrique
3.2.1 Présentation du cas test mono-cylindre et de la jupe dentelée
3.2.1.1 Géométrie d’étude et conditions de simulation
3.2.1.2 Topologie de l’écoulement de la configuration non contrôlée
3.2.2 Application de la stratégie zonale de modélisation de paroi par forçage direct
3.2.2.1 Description des étapes de la mise en oeuvre numérique
3.2.2.2 Évaluation de la précision de la formulation du terme source
3.2.3 Évaluation de la stratégie zonale de modélisation de paroi par forçage direct
3.2.3.1 Études des champs instantané et statistique
3.2.3.2 Analyse de l’impact de la jupe sur la dynamique de l’écoulement
3.3 Synthèse
4 Application de la stratégie couplée ZIBC/ZDES pour la prévision des charges latérales
4.1 Présentation du cas test et mise en oeuvre de la stratégie ZIBC
4.1.1 Présentation de la géométrie d’étude, des conditions d’essai et du maillage
4.1.2 Mise en oeuvre de l’approche numérique couplée ZIBC
4.2 Analyse de l’effet de la modélisation IBC sur les champs instantané et statistique
4.2.1 Étude de l’impact du terme source sur le champ instantané
4.2.2 Analyse de l’écoulement moyen dans le champ et à la paroi
4.3 Examen du champ fluctuant à la paroi de l’émergence
4.3.1 Analyse des signaux temporels de pression fluctuante
4.3.2 Reconstruction des efforts latéraux
4.4 Discussions sur la prévision des charges latérales en présence de frontières immergées
4.4.1 Influence de la méthode de modélisation de paroi par forçage sur la prévision des fluctuations de pression
4.4.2 Considérations sur les méthodes de reconstruction des efforts s’exerçant sur une paroi immergée
4.5 Synthèse
Conclusions
Télécharger le rapport complet
