r7.01.34 Schémas volumes finis SUSHI pour la modélisation des écoulements insaturés miscibles#
Résumé:
Cette note présente brièvement les schémas volumes finis développés dans Code_Aster. Un bref rappel des équations de comportement concernées est effectué puis les schémas utilisés sont présentés.
Table des matières
Présentation du problème : Hypothèses, Notations#
L’ensemble du modèle physique hydraulique est présenté en détail dans R7.01.11. On se contente ici de rappeler les principales variables et hypothèses ainsi que les équations traitées. Le formalisme est essentiellement issu des travaux de Coussy [5] . Les routines gérant les lois de comportement ne sont pas impactées par l’utilisation du schéma volume fini.
Cadre de modélisation#
Ces développements ne concernent que les modèles à 2 phases (liquide gaz) et 2 composants miscibles (par exemple eau et hydrogène). Nous nous plaçons donc dans le cadre des modélisations *_HH2 (cf. U2.04.05). Plus précisément on parlera ici de modélisation HH2SUDA (cf. section 4 ).
On rappelle que les 2 lois de comportement hydrauliques que l’on peut utiliser sont alors :
LIQU_AD_GAZ_VAPE : 2 composants par phase
LIQU_AD_GAZ : 2 composants dans la phase liquide, un seul dans la phase gaz (vapeur négligée)
Dans la suite on ne parlera que du modèle complet à 2 composants, 2 phases. Le cas sans vapeur revient juste à annuler les termes la concernant.
Notations#
Nous supposons que les pores du solide sont occupés par deux constituants notés \(w\) (pour l’eau) et \(h\) (pour l’hydrogène), chacun coexistant dans deux phases au maximum, l’une liquide notée \(l\) et l’autre gazeuse notée \(g\) .
Les grandeurs \(A\) associées à la phase \(p\) (\(p=l,g\) ) du composant \(c\) seront notées: \({X}_{p}^{c}\) . Concrètement, cela donne :
\({A}_{l}^{w}\) : grandeur \(A\) pour l’eau liquide,
\({A}_{g}^{w}\) : grandeur \(A\) pour la vapeur d’eau,
\({A}_{l}^{h}\) : grandeur \(A\) pour le composant h dissous dans le liquide,
\({A}_{g}^{h}\) : grandeur \(A\) pour le composant h sous forme gazeuse (par ex. hydrogène sec).
Les hypothèses générales effectuées sont les suivantes :
comportement anisotrope,
les gaz sont des gaz parfaits,
mélange idéal de gaz parfaits (pression totale = somme des pressions partielles),
équilibre thermodynamique entre les phases d’un même constituant.
Les différentes notations sont explicitées ci-après.
Variables d’état#
Les variables sont:
les pressions de chaque constituants \({P}_{p}^{c}\)
la température du milieu \(T\) .
Ces différentes variables ne sont pas totalement indépendantes. En effet, si l’on considère un seul constituant, l’équilibre thermodynamique entre ses phases impose une relation entre la pression de la vapeur et la pression du liquide de ce constituant. Finalement, il n’y a qu’une seule pression indépendante par constituant, de même qu’il n’y a qu’une seule équation de conservation de la masse. Le nombre de pressions indépendantes est donc égal au nombre de constituants indépendants. Le choix de ces pressions est libre (combinaisons des pressions des constituants) à condition que les pressions choisies, associées à la température, forment un système de variables indépendantes.
Nous avons choisi - et afin d’être homogène avec la formulation éléments finis - comme variables indépendantes et descriptives du milieu :
la pression totale du gaz \({P}_{\text{g}}={P}_{g}^{w}+{P}_{g}^{h}\) , (loi de Dalton)
La pression totale de liquide \({P}_{l}={P}_{l}^{w}+{P}_{l}^{h}\)
la pression capillaire \({P}_{c}={P}_{\text{g}}-{P}_{l}={P}_{g}^{w}+{P}_{g}^{h}-({P}_{l}^{w}+{P}_{l}^{h})\) .
Grandeurs caractéristiques de la phase solide#
On note :
La porosité : \(\phi\) ,
Le tenseur de perméabilité intrinsèque : \(k\)
Grandeurs caractéristiques des fluides#
On note :
La masse volumique de la phase \(p\) : \({\rho}_{\text{p}}\) , \({\rho}_{\text{p}}={\rho}_{p}^{w}+{\rho}_{p}^{h}\)
La viscosité de la phase \(p\) : \({\mu}_{\text{p}}\)
La saturation de la phase \(p\) : , \({S}_{l}+{S}_{\text{g}}=1\) qui est une fonction décroissante de la pression capillaire. On a donc \({S}_{l}=f({P}_{c})\) .
La perméabilité relative de la phase \(p\) : \({k}_{\text{rp}}\) fonction de la saturation
La conductivité hydraulique de la phase \(p\) : \({\lambda}_{\text{p}}^{H}\) telle que : \({\lambda}_{\text{p}}^{H}=\frac{k{k}_{\mathrm{rp}}}{{\mu}_{p}}\)
La mobilité du composant \(c,c=(h,w)\) associé à la phase \(p\) \({k}_{\text{p}}^{c}=\frac{{\rho}_{\text{p}}^{c}{k}_{\mathit{rp}}}{{\mu}_{p}}\)
La masse molaire du composant \(c\) : \({M}^{c}\)
La concentration molaire \({c}_{p}^{c}=\frac{{\rho}_{p}^{c}}{{M}^{c}}\)
La fraction massique de la phase \(p\) et du composant \(c\) : \({\zeta}_{p}^{c}=\frac{{\rho}_{p}^{c}}{{\rho}_{p}}\)
La fraction molaire : \({X}_{p}^{c}=\frac{{c}_{p}^{c}}{{c}_{p}}\) où \({c}_{p}={c}_{p}^{h}+{c}_{p}^{w}\) (dans la littérature, ces concentrations sont parfois également notées \({C}_{p}^{c}\) )
Le module de compressibilité de l’eau \({K}_{w}\)
Équations constitutives du modèle#
On ne donnera pas ici les détails permettant d’arriver aux équations finales du modèle. Pour les étapes intermédiaires, on se réfèrera à [R7.01.11]. On se contente donc ici d’un rappel succinct des principales équations.
Les équations d’équilibre sont ici données par la conservation de la masse de chaque constituant, soit :
\(\lbrace \begin{array}{}\dot{{m}_{l}^{w}}+{\dot{m}}_{g}^{w}+\text{Div}({F}_{l}^{w}+{F}_{g}^{w})=0\\ \dot{{m}_{l}^{h}}+{\dot{m}}_{g}^{h}+\text{Div}({F}_{l}^{h}+{F}_{g}^{h})=0\end{array}\)
avec \({m}_{p}^{c}\) l’apport massique du constituant \(c\) en phase \(p\) , telle que \({m}_{p}^{c}=\phi .{S}_{p}.{\rho}_{p}^{c}\) et \({F}_{p}^{c}\) le flux massique de la phase \(p\) pour le constituant \(c\) .
\({F}_{p}^{c}\) est constitué du flux massique Fickien \({J}_{p}^{c}\) et du flux massique Darcéen \({F}_{p}\) tel que :
\({F}_{p}^{c}={J}_{p}^{c}+{\rho}_{p}^{c}\frac{{F}_{p}}{{\rho}_{l}},c=(h,w);p=(l,g)\)
Lois de comportement des fluides#
Évolution de la porosité :
En l’absence de mécanique, on permet néanmoins une évolution de la porosité via le coefficient d’emmagasinement \({E}_{m}\) , telle que :
\(d\phi ={E}_{m}d{P}_{l}\)
Comportement du liquide :
On considère que l’eau peut être compressible : \(\frac{d{\rho}_{l}^{w}}{{\rho}_{l}^{w}}=\frac{{\text{dP}}_{l}^{w}}{{K}_{w}}\)
Comportement du gaz :
On considère que le gaz est soumis à la loi des gaz parfaits :
\(\frac{{P}_{g}^{c}}{{\rho}_{g}^{c}}=\frac{RT}{{M}^{c}};c=(w,h)\) où \(R\) est la constante des gaz parfaits.
Loi d’équilibre eau vapeur :
l’équilibre eau vapeur s’écrit par égalité des enthalpies libres, ce qui, pour un problème isotherme donne (confer R7.01.11) :
\(\frac{d{P}_{g}^{w}}{{\rho}_{g}^{w}}=\frac{{\text{dP}}_{l}^{w}}{{\rho}_{l}^{w}}\)
Loi d’équilibre gaz sec /dissous :
La loi de Henry relie le composant \(c\) sous sa forme gazeuse à sa forme liquide telle que :
\({P}_{g}^{h}={K}_{h}\frac{{\rho}_{l}^{h}}{{M}^{h}}\) où \({K}_{h}\) coefficient de Henry.
Remarque : on trouve souvent ce coefficient dans la littérature sous la forme \(H=\frac{1}{{K}_{h}}\)
Équations de diffusion#
Loi de Darcy :
La loi de Darcy relie le flux \({F}_{p}\) de la phase \(p\) à son gradient de pression, tel que :
\(\frac{{F}_{\text{p}}}{{\rho}_{p}}=\frac{-k{k}_{\mathrm{rp}}}{{\mu}_{p}}(\nabla {P}_{p}-{\rho}_{p}g)\)
avec \(g\) la gravité
Loi de Fick :
On écrit les flux massiques Fickiens tels que :
\({J}_{p}^{c}=-\phi {M}^{c}{S}_{p}{D}_{p}{c}_{p}\nabla {X}_{p}^{c}\) où \({D}_{p}\) est le coefficient de diffusion (où coefficient de Fick) de la phase \(p\) .
On néglige le flux Fickien de l’eau dans le liquide (concentration de l’eau dans le liquide assimilable à 1). Au final, on obtient donc :
\({F}_{l}^{w}=-{\rho}_{l}^{w}{\lambda}_{l}^{H}(\nabla {P}_{l}-{\rho}_{l}g)\)
\({F}_{l}^{h}=-{\rho}_{l}^{h}{\lambda}_{l}^{H}(\nabla {P}_{l}-{\rho}_{l}g)-\phi {M}^{h}{S}_{l}{c}_{l}{D}_{l}\nabla {X}_{l}^{h}\)
\({F}_{g}^{w}=-{\rho}_{g}^{w}{\lambda}_{\text{g}}^{H}(\nabla {P}_{g}-{\rho}_{g}g)-\phi {M}^{w}{S}_{g}{c}_{g}{D}_{g}\nabla {X}_{g}^{w}\)
\({F}_{g}^{h}=-{\rho}_{g}^{h}{\lambda}_{\text{g}}^{H}(\nabla {P}_{g}-{\rho}_{g}g)-\phi {M}^{h}{S}_{g}{c}_{g}{D}_{g}\nabla {X}_{g}^{h}\)
Les schémas volumes finis mis en œuvre dans Code_Aster#
Nous présentons dans ce chapitre les schémas volumes finis qui ont été mis en œuvre dans Code_Aster.
Généralités sur les volumes finis#
La méthode des volumes finis consiste à intégrer une ou plusieurs équations sur un volume de contrôle puis à discrétiser les flux sur chaque arête de ces volumes. Contrairement aux éléments finis, il n’y a pas a proprement parlé de formule variationnelle.
Prenons l’exemple d’une équation de la forme \(div(\Theta (X,u,\nabla u))=f(x)\) . Son approximation sur un ensemble \(\Omega \subset {ℝ}^{n}\) est une fonction constante par morceaux sur un découpage de \(\Omega\) en volumes de contrôles notés par la suite \(K\) .
Le principe de la méthode des volumes finis est d’intégrer l’équation sur chacun de ces volumes de contrôle. En utilisant la formule de Green, les équations s’écrivent :
\({\int}_{\partial K}\Theta (X,u,\nabla u).{n}_{K}d\gamma ={\int}_{K}f(x)\mathrm{dx},\forall K\)
avec \({n}_{K}\) la normale unitaire sortante du volume de contrôle \(K\) .
On voit alors qu’appliquer un schéma volumes finis consiste à approcher des flux sur des bords.
Pour une équation transitoire du type :
\(\frac{\partial u}{\partial t}+div(\Theta (X,u,\nabla u))=f(x)\) , on aura logiquement une formulation du type
\({\int}_{K}u\mathrm{dx}+{\int}_{\partial K}\Theta (X,u,\nabla u).{n}_{K}d\gamma ={\int}_{K}f(x)\mathrm{dx},\forall K\)
En résumé, une schéma volumes finis est défini par :
Le choix du volume de contrôle \(K\) qui peut être identique à la maille, ou bien une construction (de type Voronoï par exemple). La solution approchée donnée par les schémas de type volumes finis sera constante par morceaux sur ce découpage.
La position des nœuds d’approximation des inconnues pour chaque volume.
La manière d’approcher le flux sur une arête du volume de contrôle : Ils sont approchés par la formulation intégrale des équations et l’application de la formule de Green sur le bord du volume de contrôle. Lors de l’approximation, les flux doivent être conservatifs et consistants. La conservativité se traduit par la continuité des flux à l’interface, par définition assurée par les volumes finis (le flux est écrit explicitement). La consistance est atteinte lorsque la différence entre le flux réel et son approximation tend vers 0. Il existe un très grand nombre de schémas suivant la façon de calculer ce flux et la consistance requise.
L’écriture du flux peut donc relier des nœuds qui n’appartiennent pas forcement au même élément. Ils peuvent être reliés par exemple par une interface commune. On voit donc que dans ce cas, deux nœuds sont reliés, non plus quand ils appartiennent au même élément comme en éléments finis, mais quand ils appartiennent à deux voisins (reliés par une face commune).
Pour mettre en œuvre la plupart de ces schémas il faut donc avoir recours à la notion topologique de voisinage ce qui constitue une différence importante par rapport aux éléments finis.
Les volumes finis SUSHI appliqués aux écoulements diphasiques miscibles#
Les écoulements diphasiques en milieu poreux sont régis par le système d’équations :
On rappelle par ailleurs que l’expression des apports massiques est :
les flux sont alors :
Le choix habituel d’inconnues \(({P}_{c},{P}_{g})\) est retenu dans la maille ainsi que sur les arêtes de celle-ci. Ce choix est par ailleurs compatible avec le traitement unifié des écoulements miscibles et immiscibles en milieu saturé et insaturé (voir [1] ).
Le volume de contrôle correspond ici au maillage : un élément est égal à un volume de contrôle.
Les notations géométriques utilisées par la suite pour un élément \(K\) (ici triangulaire, mais le principe est le même pour des rectangle ou en 3D) sont indiquées .
On notera par la suite \(L\) le triangle voisin à \(K\) séparé par l’interface \(\sigma\) . On rappelle que la notation \(∣K∣\) désigne sa mesure (ou surface).
Illustration 1: Représentation des différentes quantités
Nous pouvons distinguer différents types de flux que nous nommerons par la suite flux «massiques» ou flux «volumiques». Les flux massiques sont en réalité les flux «physiques» classiques autrement dit ce sont les quantités \({\int}_{\partial K}-k{k}_{p}^{c}\nabla [{P}_{p}-{\rho}_{p}\mathrm{g.}x].{n}_{K,\sigma }d\gamma\) pour les termes Darcéens et \({\int}_{\partial K}-\phi {M}^{c}{S}_{p}{D}_{p}({S}_{p}){c}_{p}\nabla {X}_{p}^{c}.{n}_{K,\sigma }d\gamma\) pour les termes Fickiens.
En revanche, les flux volumiques sont pour les termes Darcéens les quantités \({\int}_{\partial K}-k\nabla [{P}_{p}-{\rho}_{p}\mathrm{g.}x],{n}_{K,\sigma }d\gamma\) et pour les termes Fickiens les quantités \({\int}_{\partial K}-\nabla {X}_{p}^{c}.{n}_{K,\sigma }d\gamma\) . Concrètement, on extrait les termes dépendant de l’espace (via la saturation) du calcul du flux volumique.
Un des choix réalisés pour l’extension du schéma SUSHI au cas non linéaire diphasique est d’approximer par des flux discrets les flux volumiques plutôt que les flux massiques. Ce choix a l’avantage de pouvoir conduire à une équation de continuité des flux volumiques linéaires. Cependant, il faut alors assurer la continuité des autres quantités discontinues pas d’autres techniques (par exemple, par décentrage ou par moyenne). Ces différentes approches sont comparées dans [1] .
On note alors la discrétisation des flux volumique telle que :
\({\int}_{\partial K}-k\nabla [{P}_{p}-{\rho}_{p}\mathrm{g.}x],{n}_{K,\sigma }d\gamma\) est discrétisé par \({\Sigma}_{\sigma \in {\epsilon}_{K}}{F}_{p,K,\sigma }({P}_{p}-{\rho}_{p}\mathrm{g.}x),p\in (l,g)\) (on somme les flux pour chaque bord \(\sigma\) appartenant à l’ensemble des arrêtes de K : \({\epsilon}_{K}\) ).
\({\int}_{\partial K}-\nabla {X}_{p}^{c}.{n}_{K,\sigma }d\gamma\) est discrétisé par \({\Sigma}_{\sigma \in {\epsilon}_{K}}{\tilde{F}}_{p,K,\sigma }({X}_{p}^{c}),c\in (h,w),p\in (l,g)\) .
Les flux volumiques discrétisés s’expriment donc de la manière suivante :
Avec
Où :
Et :
Le choix de discrétiser par la méthode des Volumes Finis SUSHI uniquement les flux volumiques conduit à devoir traiter de manière spécifique les mobilités ainsi que les termes de diffusion Fickiens. En effet, les perméabilités relatives ainsi que les tenseurs de diffusion, présents respectivement dans les mobilités et dans les termes de diffusion Fickiens, sont des quantités qui dépendent de la saturation et qui peuvent être hétérogènes.
Nous allons présenter différentes formulations possibles de la discrétisation des équations qui régissent les écoulements en milieu poreux, résultant de différents décentrages des quantités non linéaires.
On notera par la suite :
L’équation discrétisée de conservation de l’eau sur l’élément \(K\) : \({R}_{K}^{w}\)
L’équation discrétisée de conservation du composant gaz sur l’élément \(K\) : \({R}_{K}^{h}\)
L’équation discrétisée de continuité du flux d’eau à travers l’interface \(\sigma\) : \({R}_{\sigma}^{w}\)
L’équation discrétisée de continuité du flux du composant gazeux à travers l’interface \(\sigma\) : \({R}_{\sigma}^{h}\)
Le phénomène dominant dans les écoulements diphasiques en milieu poreux est la diffusion Darcéenne, nous nous focaliserons donc sur l’étude de la continuité des mobilités \({\rho}_{p}\frac{{k}_{r,p}}{{\mu}_{p}}\) . Les termes de diffusion Fickiens seront regardés en second lieu afin uniquement d’assurer leur continuité une fois celle des mobilités assurée.
Lorsque nous sommes en présence de termes non linéaires qui dépendent de l’espace, ici les mobilités, nous devons assurer leur continuité en cas de problème hétérogène. Pour cela nous pouvons les décentrer. Le décentrage des mobilités a pour objectif d’assurer la monotonie et la stabilité des schémas. Cependant, il diminue la précision des calculs en ajoutant de la diffusion numérique. Il existe de nombreux décentrages dans la littérature, le terme non linéaire de mobilité \({\rho}_{p}\frac{{k}_{r,p}}{{\mu}_{p}}\) dépendant de l’espace, pouvant être calculé par différents schémas :
le schéma amont qui décentre ce terme en fonction du sens de l’écoulement de la phase dans la maille courante ou dans la maille voisine par l’arête considérée,
le schéma causal qui va décentrer différemment ce terme lorsqu’il est facteur d’un terme de pression ou facteur d’un terme de gravité. Le décentrage se fera tout de même dans la maille courante ou dans la maille voisine,
le schéma à nombre de Péclet variable qui permet de décomposer la partie convective du flux en une combinaison linéaire entre le flux décentré du schéma amont et le flux centré. Ainsi le paramètre de cette combinaison peut être ajusté localement sur chaque arête en fonction de la diffusion introduite par la pression capillaire.
Nous présenterons ci-dessous la formulations retenue dans Code_Aster.
Initialement 3 formulations ont été implémentées dans Code_Aster: deux d’entre elles vont décentrer par le schéma amont les mobilités, soit sur l’arête de la maille courante soit sur la maille voisine par l’arête et la troisième formulation va utiliser des mobilités centrées dans la maille courante.
Nous avons décidé de conserver une unique formulation (pour des raison de stabilité et de temps de calcul): la formulation décentré arête.
Formulation : schéma Volumes Finis Décentré Arête (VFDA)
L’intérêt de cette formulation est à la fois de garder la notion de décentrage, tout en gardant un stencil de la matrice jacobienne petit. Pour cela, nous continuons d’écrire la continuité des flux totaux massiques mais nous décentrons les mobilités sur les arêtes de la maille courante. Nous considérons donc la continuité des quantités suivantes :
\({\int}_{\partial K}\left[-k{k}_{l}^{w}\nabla [{P}_{l}-{\rho}_{l}\mathrm{g.}x]-k{k}_{g}^{w}\nabla [{P}_{g}-{\rho}_{g}\mathrm{g.}x]-\phi {M}^{w}{S}_{g}{D}_{g}{c}_{g}\nabla {X}_{g}^{w}\right].{n}_{K,\sigma }d\gamma\)
et \({\int}_{\partial K}\left[-k{k}_{l}^{h}\nabla [{P}_{l}-{\rho}_{l}\mathrm{g.}x]-k{k}_{g}^{h}\nabla [{P}_{g}-{\rho}_{g}\mathrm{g.}x]-\phi {M}^{h}{S}_{g}{D}_{g}{c}_{g}\nabla {X}_{g}^{h}-\phi {M}^{h}{S}_{l}{D}_{l}{c}_{l}\nabla {X}_{l}^{h}\right].{n}_{K,\sigma }d\gamma\)
Le décentrage s’écrit alors :
En ce qui concerne les termes de diffusion Fickiens, étant donné que la diffusion moléculaire n’est pas le phénomène dominant, nous décidons de ne pas leur appliquer de traitement particulier et nous prendrons donc leurs valeurs dans la maille courante.
Ainsi nous pourrons, dans cette formulation, ne pas faire la distinction entre les arêtes internes et les arêtes de bords. Un bémol à cette formulation pourrait être qu’au final nous assurons doublement la continuité des flux massiques Darcéens puisqu’ils sont continus grâce aux équations de la continuité des flux mais également grâce à leur écriture dans la conservation de la masse.
Malgré cela cette formulation semble pouvoir être performante. En effet, elle n’utilise pas les voisins des arêtes, ce qui lui permet d’avoir un temps de résolution faible. Par ailleurs, elle permet de résoudre des problèmes hétérogènes.
La discrétisation des deux premières équations du système s’écrit de la manière suivante :
Mise en œuvre dans Code_Aster#
Dans ce chapitre, nous précisons comment sont intégrées les relations décrites au chapitre 3.
A ttention les volumes finis ne sont actuellement disponibles qu’avec la loi de couplage hydraulique HYDR_VGMqui permet de gérer correctement l’apparition/disparation de phase en traitant des pressions capillaires négatives. Les seuls modèles hydrauliques aujourd’hui disponibles sont donc du type Mualem Van-Genuchten.
Description des éléments#
2 types de modélisations existent suivant le schéma utilisé :
Modélisation |
Schéma correspondant |
Loi d’écoulement compatible |
D_PLAN_HH2SUDA |
VFDA |
LIQU_AD_GAZ_VAPE, LIQU_AD_GAZ |
3D_HH2SUDA |
VFDA |
LIQU_AD_GAZ_VAPE, LIQU_AD_GAZ |
Tableau 4.1-1 : Modèles volumes finis
Les inconnues principales sont la pression capillaire (PRE1) et la pression de gaz (PRE2) et se situent au centre des mailles ainsi qu’aux milieux des arêtes (cf. ).
Illustration 2: Elément quadratique
Ce qui donne donc pour un quadrangle un stockage de la sorte :
Support |
DDL |
Face \(\sigma 1\) |
PRE1 |
PRE2 |
|
Face \(\sigma 2\) |
PRE1 |
PRE2 |
|
Face \(\sigma 3\) |
PRE1 |
PRE2 |
|
Face \(\sigma 4\) |
PRE1 |
PRE2 |
|
Centre \(K\) |
PRE1 |
PRE2 |
Tableau 4.1-2 : Stockage des inconnues
Les éléments mailles sont définis en 2D pour des triangles à 7 nœuds et des quadrilatères à 9 nœuds (sommets «inutilisés»+ milieux des arrêtes + centre) ainsi qu’en 3D pour des hexaèdres à 27 nœuds.
On ne dispose pas d’éléments de maillages qui excluraient les nœuds sommets mais ces derniers ne sont pas pris en compte (cf. ).
Calcul des contraintes et déformations généralisées#
Les schémas volumes finis bénéficient largement de la structure définie pour les éléments finis en [R7.01.10] et [R7.01.11].
Ainsi les contraintes généralisées au centre de l’élément sont physiquement les mêmes qu’en éléments finis, à savoir :
\({m}_{l}^{w},{\mathrm{F}}_{l}^{\mathrm{w}};{m}_{g}^{w},{\mathrm{F}}_{\mathrm{g}}^{\mathrm{w}};{m}_{g}^{h},{\mathrm{F}}_{\mathrm{g}}^{\mathrm{h}};{m}_{l}^{h},{\mathrm{F}}_{l}^{\mathrm{h}}\) ainsi que les déformations généralisées: \({p}_{c},\nabla {p}_{c};{p}_{\text{g}},\nabla {p}_{\text{g}}\) .
Aux interfaces en revanche, les flux contiennent en fait le nécessaire à l’équation de continuité. Ce dernier peut être différent selon le schéma utilisé (voir section 3.2.4 ).
Nom de composante Aster |
Contenu au centre de \(K\) |
Contenu sur les faces \(\sigma \epsilon \delta K\) |
M11 |
\({\rho}_{l}^{w}\varphi {S}_{l}-{({\rho}_{l}^{w}\varphi {S}_{l})}^{-}\) |
0 |
FH11* |
\(\sum_{\delta K}^{}{F}_{l}^{w}.n\) |
\({F}_{l}^{w}.{n}_{K,\sigma }+{F}_{g}^{w}.{n}_{K,\sigma }\) |
M12 |
\({\rho}_{g}^{w}\varphi {S}_{g}-{({\rho}_{g}^{w}\varphi {S}_{g})}^{-}\) |
0 |
FH12* |
\(\sum_{\delta K}^{}{F}_{g}^{w}.n\) |
\({F}_{l}^{h}.{n}_{K,\sigma }+{F}_{g}^{h}.{n}_{K,\sigma }\) |
M21 |
\({\rho}_{g}^{h}\varphi {S}_{l}-{({\rho}_{g}^{h}\varphi {S}_{g})}^{-}\) |
0 |
FH21* |
\(\sum_{\delta K}^{}{F}_{g}^{h}.n\) |
0 |
M22 |
\({\rho}_{l}^{h}\varphi {S}_{l}-{({\rho}_{l}^{h}\varphi {S}_{l})}^{-}\) |
0 |
FH22* |
\(\sum_{\delta K}^{}{F}_{l}^{h}.n\) |
0 |
Tableau 4.2-1 : Contraintes généralisées pour les schémas VFDA (* HH2SUDA )
Intégration#
Comme en éléments finis, la boucle principale d’intégration se fait par élément. En revanche au sein d’un élément on bouclera sur les nœuds (il n’y a plus de notion de points d’intégration, mais juste des points d’approximation qui sont ici les nœuds). Ce schéma ayant ses nœuds au centre ainsi que sur chaque arête, cela permet ici implicitement de faire une boucle sur les interfaces. Finalement, la structure éléments finis n’est pas modifiée.
Au final on peut résumer l’intégration des différentes équations traitées dans le tableau suivant (écrit pour un quadrangle mais aisément généralisable) :
Support |
Équation |
Type |
Face \(\sigma 1\) |
\({R}_{\sigma 1}^{w}\) |
Continuité flux d’eau |
\({R}_{\sigma 1}^{h}\) |
Continuité flux de \(h\) |
|
Face \(\sigma 2\) |
\({R}_{\sigma 2}^{w}\) |
Continuité flux d’eau |
\({R}_{\sigma 2}^{h}\) |
Continuité flux de \(h\) |
|
Face \(\sigma 3\) |
\({R}_{\sigma 3}^{w}\) |
Continuité flux d’eau |
\({R}_{\sigma 3}^{h}\) |
Continuité flux de \(h\) |
|
Face \(\sigma 4\) |
\({R}_{\sigma 4}^{w}\) |
Continuité flux d’eau |
\({R}_{\sigma 4}^{h}\) |
Continuité flux de \(h\) |
|
Centre \(K\) |
\({R}_{K}^{w}\) |
Conservation de la masse de l’eau |
\({R}_{K}^{h}\) |
Conservation de la masse de \(h\) |
Variables internes#
Les variables internessont ici :
Numéro |
Nom de composante Aster |
Contenu |
1 |
\(\mathit{V1}\) |
\({\rho}_{\text{lq}}-{\rho}_{{\text{lq}}_{}}^{0}\) |
2 |
\(\mathit{V2}\) |
\(\phi -{\phi }^{0}\) |
3 |
\(\mathit{V3}\) |
\({p}_{\text{vp}}-{p}_{\text{vp}}^{0}\) |
4 |
\(\mathit{V4}\) |
\({S}_{\text{lq}}\) |
5 |
\(\mathit{V5}\) |
\({P}_{\text{c}}\) |
6 |
\(\mathit{V6}\) |
\({P}_{\text{g}}\) |
Remarque : on rappelle qu’ici les points de Gauss au sens Aster sont les nœuds de l’élément.
Validation#
Le tableau suivant présente quelques exemples de cas tests de validation pour des problèmes physiques classiques :
Cas test |
Phénomène |
Modélisations testées |
wtnp117 |
Rééquilibrage capillaire |
D_PLAN_HH2SUDA(c) |
wtnp120 |
Apparition/disparition de phase dans un barreau |
D_PLAN_HH2SUDA(a) 3D_HH2SUDA(c) |
wtnp121 |
Barreau saturé en eau soumis à un choc de pression |
D_PLAN_HH2SUDA(a) 3D_HH2SUDA(e) |
wtnp122 |
Barreau saturé en gaz soumis à un choc de pression |
D_PLAN_HH2SUDA(a) |
wtnp123 |
Injection de gaz autours d’une galerie |
D_PLAN_HH2SUDA(a) |
wtnp124 |
Test de Liakopoulos : drainage gravitaire d’une colonne d’eau |
D_PLAN_HH2SUDA(a) |
Bibliographie#
Angelini : «Étude des schémas numériques pour les écoulements diphasiques en milieux poreux déformables pour des maillages quelconques. Application au stockage des déchets radioactifs. Thèse d’Ophélie Angélini. » Note H-T64-2011-01146
Granet : «Chantier numérique en THM dans Code_Aster» Note HT-64-2005-01621
Eymard R., Gallouët, T., Herbin R. (2000) «Finite Volume Method», Handbook of Numerical Analysis, vol. VII, p. 713-1020. Editors: P.G. Ciarlet and J.L. Lions.
Eymard R., Gallouët, T., Herbin R. (2007) «A New Finite Volume Scheme fpr anisotropic diffusion problems on general grids : convergence analysis», C.R., Math, Acoad. Sci. Paris, vol 344, no. 6, pp.403-406.
Coussy: «Mécanique des milieux poreux». éditions TECHNIP.
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
02/10/10 |
S.Granet EDF-R&D/AMA |
Version initiale |