Décomposition de Hodge-Helmholtz Discrète

21ème Congrès Français de Mécanique Bordeaux, 26 au 30 août 2013 Décomposition de Hodge-Helmholtz discrète A. Lemoinea , J.-P. Caltagironea , M. Azaïeza , S. Vincenta a. Univ. Bordeaux, I2M, UMR 5295, F-33400 Talence, France. Résumé : La quasi omniprésence des champs de vecteurs dans les divers domaines de la physique en fait un objet d’étude à part entière. La présence de structures imbriquées (tourbillons, sources, puits) au sein de ces champs de vecteurs peut en rendre l’analyse difficile. Dans cet article, nous nous intéressons à une décomposition particulière de ces champs permettant d’en faire ressortir les diverses structures : la décomposition de Hodge-Helmholtz. Cette décomposition se fera dans un cadre numérique, donc discret. Nous présentons dans cet article l’établissement d’un outil polyvalent de décomposition discrète, utilisant une méthode de type mimétique pour la discrétisation. La finalité étant d’évaluer les capacités de la méthode à analyser la turbulence et à améliorer les algorithmes numériques de résolution des équations de Navier-Stokes. Abstract : Vector fields are present in almost every branch of physics. The presence of nested structures, such as vortices, sources or sinks, can make analysis difficult. We present here a decomposition of vector fields able to untangle these structures : the Helmholtz-Hodge decomposition. This decomposition will take place in a discrete framework. We present here a versatile tool able to manage a discrete HelmholtzHodge decomposition using a mimetic method. We aim to evaluate the abilities of this decomposition to analyse the turbulence and we aim to improve some numerical algorithms for solving the Navier-Stokes equations. Mots clefs : 1 Hodge-Helmholtz, Turbulence, Mimetic Introduction La décomposition de Hodge-Helmholtz (DHH ) est une décomposition de champs de vecteurs en la somme de trois composantes (1) : une composante irrotationnelle uϕ issue du gradient d’un potentiel scalaire ϕ, une composante solénoïdale issue du rotationnel d’un potentiel vecteur ψ et une composante ne pouvant s’intégrer dans l’une des deux composantes précédentes, à la fois irrotationnelle et solénoïdale, dite harmonique uh . u = = uϕ grad ϕ | {z } irrotationnel + + uψ rot ψ | {z } solénoïdal + + uh uh |{z} (1) harmonique Les champs de vecteurs apparaissent dans de nombreuses branches de la physique. La complexité de ces champs peut rendre leur étude difficile. Par exemple, en mécanique des fluides, les écoulements turbulents présentent des motifs très complexes, composés de structures imbriquées, telles que des tourbillons, des sources et des puits. La DHH permet de séparer ces différentes structures et ainsi d’en faciliter l’analyse. Notons d’ailleurs que cette décomposition a été originellement introduite par Helmholtz [4] pour étudier les écoulements de fluides parfaits. D’autres applications plus exotiques sont disponibles dans la littérature. On citera l’application de la DHH par Palit et al. [7] à la détection de motifs dans des images, tel que le suivi d’un ouragan à partir d’images satellite ou la détection du point de courbure maximale d’empreintes digitales. Le lecteur pourra se référer à la revue de l’état 1 21ème Congrès Français de Mécanique Bordeaux, 26 au 30 août 2013 de l’art réalisée récemment par de Bhatia et al. [2] qui recense les nombreuses applications de cette décomposition. L’objectif de cette contribution est l’évaluation des capacités de la DHH à analyser des champs turbulents discrets et à améliorer les algorithmes de résolution numérique des équations de Navier-Stokes. L’accomplissement de ces objectifs nécessite la définition et la mise en œuvre d’algorithmes polyvalents de décomposition discrète, ce qui est l’objet de cet article. Dans une première partie, nous rappellerons les propriétés de la DHH qui nous conduiront au choix de la méthode numérique qui sera présentée dans une seconde partie. Nous présenterons ensuite deux méthodes de décomposition, puis nous présenterons des résultats en mettant en avant l’importance du choix des conditions aux limites. 2 Propriétés Dans le cadre de cette étude, la décomposition sera effectuée sur un domaine de Rn (n ∈ {2, 3}) noté Ω. La DHH est une décomposition de l’espace fonctionnel (L2 (Ω))n telle que présentée dans l’équation (1). Les extrema des potentiels correspondent au centre des points critiques des composantes correspondantes. Ainsi, un extremum dans le potentiel scalaire correspondra à un puits, une source ou un point-selle dans la composante irrotationnelle. Un extremum dans le potentiel vecteur correspondra au centre d’un vortex. L’orthogonalité de cette décomposition, impliquant son unicité, est pilotée par les conditions aux limites imposées aux composantes ou aux potentiels. La multitude des conditions aux limites donne des décompositions aux propriétés physiques encore mal comprises par la communauté. Bhatia et al. [2] soulignent ce dernier point en évoquant l’existence chez les auteurs traitant du sujet une scission concernant la pertinence physique de l’orthogonalité de cette décomposition. Établir un outil capable de manipuler un maximum de conditions aux limites est primordial pour les besoins de notre étude. Nous souhaitons également que la composante irrotationnelle soit exactement à rotationnel nul et que la composante solénoïdale soit réellement à divergence nulle. Ces propriétés se traduisent de manière vectorielle par deux identités sur les opérateurs différentiels : rot(grad) = 0 et div(rot) = 0. Ces propriétés doivent également être vérifiées au niveau discret. Nombre d’auteurs (Hyman & Shashkov [5], Kreeft et al. [6], Tong et al. [8]) s’attachent à vérifier ces propriétés dans leurs méthodes numériques. La méthode numérique de discrétisation que nous avons choisi a été sélectionnée pour son respect de ces identités vectorielles : la méthode mimétique. 3 Méthode mimétique Les méthodes mimétiques sont des méthodes numériques mimant les propriétés des opérateurs différentiels au niveau discret. Ces méthodes reposent sur une représentation discrète pertinente des grandeurs physiques sur les divers éléments composant le domaine discret. Le domaine discret est un maillage polyédrique (ou polygonal en 2D) composé de points, d’arêtes, de faces et de cellules orientés, tels que présentés sur la figure 1. Il existe deux types d’orientations. La première, dite interne, indépendante de l’espace ambiant dans lequel est plongé l’élément de maillage. Ce sera ce type d’orientation qui sera utilisé pour les points et les arêtes. Les points peuvent être considérés soit comme des sources, soit comme des puits. Les arêtes peuvent être parcourues dans un sens ou dans l’autre. Le second type d’orientation est dite externe. Elle dépend de l’espace ambiant dans lequel est plongé le maillage. Ainsi, les faces seront orientées par leurs normales et les cellules seront considérées soit comme des sources, soit comme des puits. Alors que les volumes finis et les éléments finis utilisent seulement les points et les cellules pour représenter les champs discrets, les méthodes mimétiques utilisent tous les éléments du maillage. Les champs scalaires seront représentés sur les points et les cellules, alors que les champs vectoriels préféreront les arêtes et les faces. Par exemple, un potentiel sera représenté sur les points, une circulation sur les arêtes, un flux à travers une surface et une densité dans un volume. Ainsi, une grandeur physique discrète ne sera pas localisée en un point, mais sur son élément d’intégration. Le lecteur pourra trouver chez Tonti [9] la signification physique des représentations discrète des grandeurs. 2 21ème Congrès Français de Mécanique (a) Point Potentiel (b) Arête Circulation Bordeaux, 26 au 30 août 2013 (c) Face Flux (d) Cellule Densité Figure 2 – Maillage primal (en gras) et maillage dual barycentrique Figure 1 – Éléments composant le maillage primal (en pointillés) Ces méthodes nécessitent également la définition d’un maillage dual. À l’instar du maillage primal, le maillage dual est composé d’éléments orientés. Les sommets de ces éléments sont construits à partir des barycentres de chacun des éléments du primal. À chaque élément du maillage primal est associé un unique élément du maillage dual. À un point primal sera associé une cellule duale, à une arête duale sera associé une face duale, etc. La figure 2 représente un maillage bidimensionnel. Le maillage primal y est représenté en traits continus et le maillage dual en pointillés. Le maillage et les champs étant définis, nous pouvons définir les opérateurs discrets. Ces opérateurs découlent des formules de Stokes (2). P, A, F et C désignent respectivement un point, une arête, une face et une cellule. t désigne la tangente d’une arête et n le vecteur normal à une face. p est un champ scalaire, ψ et u sont des champs vectoriels. Z grad p · t dl = p(P2 ) − p(P1 ) A Z rot u · n ds = F Z u · t dl ∂F Z div ψ dv = C Z ψ · n ds (2) ∂C Grâce à ces formules, le résultat d’un opérateur discret défini sur des éléments de maillage de dimension k + 1 se résume à sommer la valeur des champs définis sur le bord de ces éléments, c’est-à-dire des éléments de dimension k. Par exemple, étant donnée une cellule, la divergence correspond à la somme des flux à travers ses faces. Cette définition des opérateurs discrets permet d’assurer naturellement le respect des identités vectorielles énoncées auparavant. Les champs étant définis sur les bons éléments, aucune approximation n’est nécessaire à l’application de ces formules. Elles sont, du fait, exactes. Nous pouvons également définir ces opérateurs sur le maillage dual. En pratique, les équations nécessitent des opérateurs différentiels d’ordre 2, comme le laplacien. Le laplacien scalaire est la composition d’un gradient par une divergence. Dans les méthodes mimétiques, l’un de ces opérateurs est défini sur le maillage dual et l’autre sur le maillage primal. L’assemblage de l’opérateur final passe par une phase de transfert d’informations entre ces deux maillages. En termes physiques, cette opération correspond à l’application d’une équation constitutive (voir Tonti [9]), comme l’égalisation d’un flux et d’une circulation, éventuellement modulée par un coefficient de diffusion. En termes mathématiques cette opération revient à appliquer un opérateur de Hodge discret. Toute la qualité de la méthode numérique, notamment l’erreur et l’ordre de convergence, repose sur la définition de ce transfert. Les conditions aux limites naturellement disponibles pour ce genre de méthode seront de même nature que les grandeurs physiques présentes sur le maillage (potentiel, circulation, flux), avec la possibilité de les appliquer sur le maillage primal ou le maillage dual. 4 Algorithmes de décomposition Nous distinguons deux types de décompositions. Celles qui extraient les potentiels et celles qui extraient les composantes. Chacune de ces configurations nécessite trois étapes pour obtenir la décomposition, une relative à chaque composante. Les parties relatives aux composantes irrotationnelles et solénoïdales nécessitent la résolution d’un système linéaire. La composante harmonique est ensuite extraite par soustraction. 3 21ème Congrès Français de Mécanique 4.1 Bordeaux, 26 au 30 août 2013 Extraction des potentiels L’extraction des potentiels est la plus répandue des méthodes. Cette méthode est décrite dans le livre de Girault & Raviart [3]. Elle est composée d’un laplacien et d’un système rot-rot (3). div(grad ϕ) = div u (3) rot(rot ψ) = rot u D’un point de vue numérique, le champ initial u peut être discrétisé soit comme un flux, soit comme une circulation. De là dépend la nature des potentiels résultants. Dans le cas où le champ initial est discrétisé comme un flux, le potentiel scalaire ϕ sera un potentiel défini sur les points et le potentiel vecteur ψ sera un flux défini sur les faces. Dans l’autre cas, le potentiel scalaire sera une densité définie sur les cellules et le potentiel vecteur sera une circulation définie sur les arêtes. Les conditions aux limites naturellement applicables au laplacien sont l’imposition du flux normal (Neumann) ou l’imposition de la valeur (Dirichlet). L’opérateur rot-rot admet une condition de type circulation ψ · t et une condition qui se peut se traduire par rot ψ · t, où t désigne le vecteur tangent aux arêtes primales ou duales. Notons que les composantes sont obtenues en appliquant le gradient et le rotationnel discret aux potentiels. 4.2 Extraction des composantes L’extraction des composantes est une méthode originale proposée par Angot et al. [1]. Elle consiste à obtenir les composantes de la décomposition sans rechercher les potentiels. Elle est composée de deux systèmes d’équations (4). La première, dénommée Vector Penalty Projection (VPP ) et la seconde, dénommée Rotational Penalty Projection (RPP ). Elles dépendent d’un paramètre de pénalisation ε qui doit être choisi le plus petit possible. εuϕ + grad(div uϕ ) = grad(div u) εuψ + rot(rot uψ ) = rot(rot u) (4) Comme pour la méthode d’extraction des potentiels, le champ de vecteurs u peut être discrétisé soit comme une circulation, soit comme un flux. Dans les deux cas, la nature discrète des composantes résultantes uϕ et uψ doit être la même que celle de u. Les conditions aux limites naturellement applicables à VPP sont le flux normal uϕ · n sur les faces primales ou duales, la divergence div uϕ |∂Ω sur les points duaux et la circulation uϕ · t sur les arêtes primales. RPP admet naturellement des conditions aux limites du type circulation uψ · t sur les arêtes primales ou duales, le flux normal uψ · t sur les faces primales et une condition du type rot uψ · t sur les arêtes duales. L’admissibilité de telle ou telle condition limite dépend de la nature discrète initiale du champ u (flux ou circulation). 5 Résultats Afin de simplifier les calculs et la lecture des résultats, nous utiliserons ici un maillage bidimensionnel cartésien. L’avantage de ce type de maillage est qu’il ne nécessite qu’un opérateur de Hodge s’écrivant sous la forme d’une matrice diagonale tout en permettant d’obtenir une erreur de convergence à l’ordre 2. De plus, nous proposons dans cette partie d’évaluer qualitativement les résultats de la décomposition de Hodge-Helmholtz discrète pour différents types de conditions aux limites. Considérons un domaine carré Ω = [−1, 1]2 ainsi qu’un champ initial u égal à la somme des deux composantes présentées dans l’équation (5). uϕ (x, y) =  sin(πx) cos(πy) cos(πx) sin(πy)  uψ (x, y) =  cos(πx) sin(πy) − sin(πx) cos(πy)  (5) Comme représenté sur la figure (3a), le champ initial se compose de tourbillons, de puits et de sources situés au centre, aux quatre coins et aux centres des côtés du domaine. 4 21ème Congrès Français de Mécanique (a) u original Bordeaux, 26 au 30 août 2013 (b) uϕ original (c) uψ original Figure 3 – Champs de vecteurs initiaux Supposons maintenant que nous ne connaissons pas les composantes initiales. Appliquons VPP et RPP au champ u avec une condition limite tangentielle nulle pour la composante uϕ et une condition limite tangentielle non homogène pour la composante uψ : uψ · t = u · t. Notons que dans ce cas, la composante harmonique est nulle. Les résultats sont présentés sur les figures 4b et 4c. Le maillage est une grille cartésienne contenant 642 points. L’inversion du système a été réalisé par une méthode BiCGStab(2). Le paramètre de pénalisation a été imposé à ε = 10−15 . Le rotationnel discret de uϕ et la divergence discrète et uψ sont nuls à l’erreur machine près. (a) u original (b) uϕ , uϕ · t = 0 sur ∂Ω (c) uψ , uψ · t = u · t sur ∂Ω Figure 4 – Influence des conditions aux limites sur la décomposition de Hodge-Helmholtz discrète. On peut noter une différence de position des puits et des sources dans la composante irrotationnelle par rapport à la position de ces derniers dans le champ original. Les conditions aux limites n’étant pas celles des composantes originales, les composantes obtenues sont différentes. D’un point de vue qualitatif, les positions des points critiques dans la composante irrotationnelle s’en trouvent altérées. En particulier, on remarque que plus on s’approche des bords, plus la composante obtenue diffère de la composante originale. Pire encore, les points critiques des composantes ne se superposent pas avec ceux du champ initial. Ce genre de condition limite se prête donc mal à la détection de structures. Dans le cas de champs périodiques, le problème des conditions aux limites ne se pose pas. Remarquons que le champ u est 2π-périodique. Sa décomposition discrète en considérant le domaine périodique nous permet de retrouver les composantes originales. 5 21ème Congrès Français de Mécanique 6 Bordeaux, 26 au 30 août 2013 Conclusion et perspectives Nous avons mis en place un outil de décomposition de Hodge-Helmholtz discrète respectant les identités vectorielles grâce à une méthode mimétique, garantissant le respect des propriétés de la décomposition au niveau discret. Deux types d’algorithmes de décomposition ont été présentés, l’un permettant d’extraire les potentiels, l’autre permettant d’extraire directement les composantes. Grâce à cet outil, nous pouvons observer le comportement de la décomposition face aux différentes conditions aux limites. Les expériences numériques ont montré que toutes les conditions aux limites ne sont pas équivalentes du point de vue des propriétés des décompositions qu’elles engendrent. Ainsi, nous avons vu un exemple de décomposition où les points critiques (sources et puits) de la composante irrotationnelle ne coïncidaient pas avec avec ceux du champ initial, concluant que les conditions aux limites utilisées n’étaient pas aptes à générer une décomposition adaptée à la détection de structures. Malgré cela, cette décomposition reste bien une décomposition de Hodge-Helmholtz. L’analyse exhaustive des diverses conditions aux limites fera l’objet d’un travail d’étude approfondi. L’objectif étant de déterminer quelles sont les conditions aux limites les plus pertinentes pour l’étude de la physique. Par ailleurs, dans le cas d’un domaine périodique, ce problème ne se présente pas. La décomposition de champs périodiques sera utile pour l’étude de champs turbulents homogènes isotropes compressibles dans le but de faire émerger de nouveaux résultats. On pourra, par exemple, envisager d’étudier l’évolution de la distribution d’énergie dans les composantes de la décomposition en fonction du temps. Références [1] Philippe Angot, Jean-Paul Caltagirone, and Pierre Fabrie. Fast discrete Helmholtz-Hodge decompositions in bounded domains. Applied Mathematics Letters, 26, 2013. 8 pages. [2] Harsh Bhatia, Gregory Norgard, Valerio Pascucci, and Peer-Timo Bremer. The Helmholtz-Hodge Decomposition - A Survey. IEEE Transactions on Visualization and Computer Graphics, 99(PrePrints) :1, 2012. [3] V. Girault and P.A. Raviart. Finite element methods for Navier-Stokes equations : theory and algorithms. Springer series in computational mathematics. Springer-Verlag, 1986. [4] Hermann L. von Helmholtz. Über integrale der hydrodynamischen gleichungen, welche den wirbelbewegungen entsprechen. Journal für die reine und angewandte Mathematik, 55 :25–55, 1858. [5] James M. Hyman and Mikhail J. Shashkov. The orthogonal decomposition theorems for mimetic finite difference methods. SIAM J. Numer. Anal., 36(3) :788–818 (electronic), 1999. [6] Jasper Kreeft, Artur Palha, and Marc Gerritsma. Mimetic framework on curvilinear quadrilaterals of arbitrary order. Submitted to Foundations of Computational Mathematics, arXiv :1111.4304, November 2011. [7] Biswaroop Palit, Anup Basu, and Mrinal K. Mandal. Applications of the Discrete Hodge-Helmholtz Decomposition to Image and Video Processing. Lecture Notes in Computer Science, 3776/2005 :497– 502, 2005. [8] Yiying Tong, Santiago Lombeyda, Anil N Hirani, and Mathieu Desbrun. Discrete multiscale vector field decomposition. In ACM SIGGRAPH 2003 Papers on - SIGGRAPH ’03, number 1, page 445, New York, New York, USA, 2003. ACM Press. [9] Enzo Tonti. The reason for analogies between physical theories. Appl. Math. Modelling, 1(1) :37–50, 1976/77. 6