u2.04.05 Notice d’utilisation du modèle THM#

Résumé:

Les modélisations THM traitent exclusivement de l’évolution des comportements Thermo-Hydro-Mécaniques des milieux poreux saturés, ou non, par un ou deux fluides. La description numérique et physique des ces modélisations dans Code_Aster est détaillée dans [R7.01.10] et [R7.01.11].

On détaille dans cette documentation la procédure à suivre pour réaliser des simulations dans le contexte des modélisations THM. On décrit dans la première partie les différentes étapes de calculs dans le cadre général des milieux continus dits «généralisés» (cf [R7.01.10]). Dans cette partie, on ne détaillera pas ce qui concerne les modèles mécaniques décrits ailleurs. On restreint dans la seconde partie l’application de ces modélisations pour traiter les milieux poreux qui subissent une dégradation de leur propriétés mécaniques par endommagement du sol ou des roches par exemple. Pour cela, on étend les modélisations THM aux milieux à microstructure en prenant en compte les effets second gradient (cf [R5.04.03]). L’objectif est de corriger la dépendance à la discrétisation spatiale des solutions lorsque la loi de comportement mécanique considérée est de type adoucissante – ce qui est le cas pour tout matériau fragile, et donc les sols en particulier.

Les différentes étapes d’un calcul THM#

Choix du modèle#

Le traitement numérique en THM nécessite un maillage quadratique puisque les éléments sont de type \(\text{P2}\) en déplacement et \(\text{P1}\) en pression et température afin d’éviter des problèmes d’oscillations.

Le choix se fait par l’utilisation de la commande AFFE_MODELE comme dans l’exemple ci-dessous:

MODELE=AFFE_MODELE(MAILLAGE=MAIL,
                   AFFE=_F(TOUT='OUI',
                           PHENOMENE='MECANIQUE',
                           MODELISATION='AXIS_THH2MD',),
                   )

Dans tous les cas, le phénomène est ``MECANIQUE``(même si la modélisation ne contient pas la mécanique).

L’utilisateur doit renseigner ensuite de manière obligatoire le mot clé MODELISATION. Ce mot clé permet de définir le type d’élément affecté à un type de maille. Les modélisations disponibles en THM sont indiqués dans le tableau suivant.

Remarque

Les modélisations se terminant par la lettre D indiquent que l’on fait un traitement permettant de diagonaliser (« lumper ») la matrice afin d’éviter les oscillations pour les problèmes hydrauliques. Pour cela les points d’intégration sont pris aux sommets des éléments. Ce traitement étant peu adapté à la mécanique, on dispose également d’une modélisation dite «sélective». Dans ce cas, les termes capacitifs sont intégrés aux sommets alors que les termes diffusifs sont intégrés aux points de Gauss. Ces modélisations se terminent par un S. Les autres modélisations intègrent tout aux point de Gauss.

On conseille vivement à l’utilisateur d’utiliser les modélisations D ou S dans les cas sans mécanique et d’utiliser la modélisation S et pour les modélisations avec mécanique.

Les modélisations « classiques » (sans D ni S) sont amenées à être résorbées et sont déconseillées.

Tableau 6 Modélisations THM#

MODELISATION

Modélisation géométrique

Phénomènes pris en compte

D_PLAN_HM

plane

Mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

D_PLAN_HMD

plane

Mécanique, hydraulique avec une pression inconnue (rigidité lumpée). Modélisation valide en dynamique

D_PLAN_HMS

plane

Mécanique, hydraulique avec une pression inconnue (sélective). Modélisation valide en dynamique

D_PLAN_HM_SI

plane

Mécanique, hydraulique avec une pression inconnue (sous-intégrée)

D_PLAN_HHM

plane

Mécanique, hydraulique avec deux pressions inconnues. Modélisation valide en dynamique

D_PLAN_HHMD

plane

Mécanique, hydraulique avec deux pressions inconnues (rigidité lumpée). Modélisation valide en dynamique

D_PLAN_HHMS

plane

Mécanique, hydraulique avec deux pressions inconnues (sélective). Modélisation valide en dynamique

D_PLAN_HH2MD

plane

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

D_PLAN_HH2MS

plane

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

D_PLAN_HH2M_SI

plane

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sous-intégrée). Modélisation valide en dynamique.

D_PLAN_THHD

plane

Thermique, hydraulique avec deux pressions inconnues (lumpé)

D_PLAN_THH2D

plane

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

D_PLAN_THH2S

plane

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

D_PLAN_THM

plane

Thermique, mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

D_PLAN_THVD

plane

Thermique, mécanique, hydraulique avec deux pressions inconnues (2 phases: eau liquide et vapeur) (lumpé)

D_PLAN_THMD

plane

Thermique, mécanique, hydraulique avec une pression inconnue (lumpé)

D_PLAN_THMS

plane

Thermique, mécanique, hydraulique avec une pression inconnue (sélective)

D_PLAN_THHMD

plane

Thermique, mécanique, hydraulique avec deux pressions inconnues (lumpé)

D_PLAN_THHMS

plane

Thermique, mécanique, hydraulique avec deux pressions inconnues (sélective)

D_PLAN_HHD

plane

hydraulique avec deux pressions inconnues (lumpé)

D_PLAN_HHS

plane

hydraulique avec deux pressions inconnues (sélective)

D_PLAN_HH2D

plane

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

D_PLAN_HH2S

plane

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

D_PLAN_THH2MD

plan

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

D_PLAN_THH2MS

plane

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

AXIS_HM

axisymétrique

Mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

AXIS_HMD

axisymétrique

Mécanique, hydraulique avec une pression inconnue (rigidité lumpée). Modélisation valide en dynamique

AXIS_HMS

axisymétrique

Mécanique, hydraulique avec une pression inconnue (sélective). Modélisation valide en dynamique

AXIS_HHM

axisymétrique

Mécanique, hydraulique avec deux pressions inconnues. Modélisation valide en dynamique

AXIS_HHMD

axisymétrique

Mécanique, hydraulique avec deux pressions inconnues (rigidité lumpée). Modélisation valide en dynamique

AXIS_HHMS

axisymétrique

Mécanique, hydraulique avec deux pressions inconnues (sélective). Modélisation valide en dynamique

AXIS_HH2MD

axisymétrique

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

AXIS_HH2MS

axisymétrique

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

AXIS_THHD

axisymétrique

Thermique, hydraulique avec deux pressions inconnues (lumpé)

AXIS_THHS

axisymétrique

Thermique, hydraulique avec deux pressions inconnues (sélective)

AXIS_THH2D

axisymétrique

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

AXIS_THH2S

axisymétrique

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

AXIS_THM

axisymétrique

Thermique, mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

AXIS_THMD

axisymétrique

Thermique, mécanique, hydraulique avec une pression inconnue (lumpé)

AXIS_THMS

axisymétrique

Thermique, mécanique, hydraulique avec une pression inconnue (sélective)

AXIS_THVD

axisymétrique

Thermique, mécanique, hydraulique avec deux pressions inconnues (2 phases: eau liquide et vapeur) (lumpé)

AXIS_THHMD

axisymétrique

Thermique, mécanique, hydraulique avec deux pressions inconnues (lumpé)

AXIS_THHMS

axisymétrique

Thermique, mécanique, hydraulique avec deux pressions inconnues (sélective)

AXIS_THH2MD

axisymétrique

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

AXIS_THH2MS

axisymétrique

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

AXIS_HHD

axisymétrique

hydraulique avec deux pressions inconnues (lumpé)

AXIS_HHS

axisymétrique

hydraulique avec deux pressions inconnues (sélective)

AXIS_HH2D

axisymétrique

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

AXIS_HH2S

axisymétrique

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

3D_HM

3D

Mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

3D_HMD

3D

Mécanique, hydraulique avec une pression inconnue (rigidité lumpée). Modélisation valide en dynamique

3D_HMS

3D

Mécanique, hydraulique avec une pression inconnue (sélective). Modélisation valide en dynamique

3D_HM_SI

3D

Mécanique, hydraulique avec une pression inconnue (sous-intégrée)

3D_HHM

3D

Mécanique, hydraulique avec deux pressions inconnues. Modélisation valide en dynamique

3D_HHMD

3D

Mécanique, hydraulique avec deux pressions inconnues (rigidité lumpée). Modélisation valide en dynamique

3D_HHMS

3D

Mécanique, hydraulique avec deux pressions inconnues (sélective). Modélisation valide en dynamique

3D_HH2MD

3D

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

3D_HH2MS

3D

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

3D_HH2M_SI

3D

Mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sous-intégrée). Modélisation valide en dynamique.

3D_THHD

3D

Thermique, hydraulique avec deux pressions inconnues (rigidité lumpée). Modélisation valide en dynamique

3D_THH2D

3D

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

3D_THH2S

3D

Thermique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

3D_THM

3D

Thermique, mécanique, hydraulique avec une pression inconnue. Modélisation valide en dynamique

3D_THMD

3D

Thermique, mécanique, hydraulique avec une pression inconnue (lumpé)

3D_THMS

3D

Thermique, mécanique, hydraulique avec une pression inconnue (sélective)

3D_THVD

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues (2 phases: eau liquide et vapeur) (lumpé)

3D_THHM

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues. Modélisation valide en dynamique

3D_THHMD

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues (lumpé)

3D_THHMS

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues (sélective)

3D_THH2MD

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

3D_THH2MS

3D

Thermique, mécanique, hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

3D_HHD

3D

hydraulique avec deux pressions inconnues (lumpé)

3D_HHS

3D

hydraulique avec deux pressions inconnues (sélective)

3D_HH2D

3D

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (lumpé)

3D_HH2S

3D

hydraulique avec deux pressions inconnues et deux constituants dans la phase gazeuse (sélective)

Les inconnues principales qui sont aussi les valeurs des degrés de liberté, sont notéesdans le cas de la modélisation la plus complète (3D thermique, mécanique, hydraulique avec deux pressions inconnues).

\({\left\lbrace u\right\rbrace }^{\mathrm{ddl}}=\left\lbrace \begin{array}{}{u}_{x}\\ {u}_{y}\\ {u}_{z}\\ {\text{PRE1}}^{\mathrm{ddl}}\\ {\text{PRE2}}^{\mathrm{ddl}}\\ {T}^{\mathrm{ddl}}\end{array}\right\rbrace\)

Le contenu de \(\mathit{PRE1}\) et \(\mathit{PRE2}\) dépend du couplage choisi et sera explicité dans la section 2.2.2 .

Suivant la modélisation choisie, seuls certains de ces degrés de liberté existent. Le tableau ci-dessous résume les degrés de liberté utilisés pour chaque modélisation.

MODELISATION

\({u}_{x}\)

\({u}_{y}\)

\({u}_{z}\)

\({\mathrm{PRE1}}^{\mathrm{ddl}}\)

\({\mathrm{PRE2}}^{\mathrm{ddl}}\)

\({T}^{\mathrm{ddl}}\)

D_PLAN_HM(_,S,D,SI)

X

X

X

D_PLAN_HHM(_,S,D)

X

X

X

X

D_PLAN_HH2M(S,D,_SI)

X

X

X

X

D_PLAN_THH(S,D)

X

X

X

D_PLAN_THH2(S,D)

X

X

X

D_PLAN_THM(_,S,D)

X

X

X

X

D_PLAN_THV(_,S,D)

X

X

X

D_PLAN_THHM(S,D)

X

X

X

X

X

D_PLAN_THH2M(S,D)

X

X

X

X

X

D_PLAN_HH(S,D)

X

X

D_PLAN_HH2(S,D)

X

X

AXIS_HM(_,S,D)

X

X

X

AXIS_HHM(_,S,D)

X

X

X

X

AXIS_HH2M(S,D)

X

X

X

X

AXIS_THH(S,D)

X

X

X

AXIS_THH2(S,D)

X

X

X

AXIS_THM(_,S,D)

X

X

X

X

AXIS_THV(_,S,D)

X

X

X

AXIS_THHM(S,D)

X

X

X

X

X

AXIS_THH2M(S,D)

X

X

X

X

X

AXIS_HH (S,D)

X

X

AXIS_HH2(S,D)

X

X

3D_HM(_,S,D,SI)

X

X

X

X

3D_HHM(_,S,D)

X

X

X

X

X

3D_HH2M(S,D,_SI)

X

X

X

X

X

3D_THH(S,D)

X

X

X

3D_THH2D

X

X

X

3D_THM(_,S,D)

X

X

X

X

X

3D_THVD(_,S,D)

X

X

X

3D_THHM(_,S,D)

X

X

X

X

X

X

3D_THH2M(S,D)

X

X

X

X

X

X

3D_HH(S,D)

X

X

3D_HH2(S,D)

X

X

Les contraintes généralisées et les variables internes sont toutes indiquées en [§Annexe 1]. Les notations utilisées sont celles définies dans [R7.01.11].

Définition du matériau#

Le matériau est défini par la commande DEFI_MATERIAU comme dans l’exemple ci‑dessous:

MATERBO=DEFI_MATERIAU(  ELAS        = _F(E=5.15000000E8, NU=0.20, ALPHA=0.,),
                        COMP_THM    = 'LIQU_AD_GAZ_VAPE',
                        THM_LIQU    = _F(
                                        RHO=1000.0,
                                        UN_SUR_K=0.,
                                        ALPHA=0.,
                                        CP=0.0,
                                        VISC=VISCOLIQ,
                                        D_VISC_TEMP=DVISCOL,
                                        ),
                        THM_GAZ     = _F(
                                        MASS_MOL=0.01,
                                        CP=0.0,
                                        VISC=VISCOGAZ,
                                        D_VISC_TEMP=ZERO,
                                        ),
                        THM_VAPE_GAZ= _F(
                                        MASS_MOL=0.01,
                                        CP=0.0,
                                        VISC=VISCOGAZ,
                                        D_VISC_TEMP=ZERO,
                                        ),
                        THM_AIR_DISS= _F(
                                        CP=0.0,
                                        COEF_HENRY=HENRY
                                        ),
                        THM_INIT    = _F(
                                        TEMP=300.0,
                                        PRE1=0.0,
                                        PRE2=1.E5,
                                        PORO=1.,
                                        PRES_VAPE=1000.0,
                                        ),
                        THM_DIFFU   = _F(
                                        R_GAZ=8.32,
                                        RHO=2200.0,
                                        CP=1000.0,
                                        BIOT_COEF=1.0,
                                        SATU_PRES=SATUBO,
                                        D_SATU_PRES=DSATBO,
                                        PESA_X=0.0,
                                        PESA_Y=0.0,
                                        PESA_Z=0.0,
                                        PERM_IN=KINTBO,
                                        PERM_LIQU=UNDEMI,
                                        D_PERM_LIQU_SATU=ZERO,
                                        PERM_GAZ=UNDEMI,
                                        D_PERM_SATU_GAZ=ZERO,
                                        D_PERM_PRES_GAZ=ZERO,
                                        FICKV_T=ZERO,
                                        FICKA_T=FICK,
                                        LAMB_T=ZERO,
                                        ),
                        )

Nous allons maintenant détailler chacun des mots clés. Nous ne nous attacherons pas ici à la partie mécanique – si mécanique il y a - qui dépend de la loi de comportement choisie. On se reportera pour cela à la documentation de DEFI_MATERIAU [U4.43.01].

Mot clé simple COMP_THM#

Permet de sélectionner dès la définition du matériau la loi de couplage THM.

GAZ: loi de comportement d’un gaz parfait c’est-à-dire vérifiant la relation \(P/\rho =\mathit{RT}/\mathit{Mv}\) où \(P\) est la pression, \(\rho\) la masse volumique, \(\text{Mv}\) la masse molaire, \(R\) la constante des gaz parfaits et \(T\) la température (cf. [R7.01.11] pour plus de détails). Pour un milieu saturé uniquement. Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous le mot clé THM_GAZ.

LIQU_SATU: loi de comportement pour un milieux poreux saturé par un seul liquide (Cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous le mot clé THM_LIQU.

LIQU_GAZ_ATM: Loi de comportement pour un milieu poreux non saturé avec un liquide et du gaz à pression atmosphérique (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés THM_LIQU et THM_GAZ.

LIQU_VAPE_GAZ: loi de comportement pour un milieu poreux non saturé eau/vapeur/air sec avec changement de phase (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés : THM_LIQU, THM_VAPE et THM_GAZ.

LIQU_AD_GAZ_VAPE: loi de comportement pour un milieu poreux non saturé eau/vapeur/air sec/air dissous avec changement de phase (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés THM_LIQU, THM_VAPE, THM_GAZ et THM_AIR_DISS.

LIQU_AD_GAZ :loi de comportement pour un milieu poreux non saturé eau/air sec/air dissous avec changement de phase (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés THM_LIQU, THM_GAZ et THM_AIR_DISS.

LIQU_VAPE: loi de comportement pour un milieux poreux saturé par un composant présent sous forme liquide ou vapeur avec changement de phase (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés THM_LIQU et THM_VAPE. **Cette loi n’est valable que pour les modélisations de type THVD .

LIQU_GAZ: loi de comportement pour un milieu poreux non saturé liquide/gaz sans changement de phase (cf. [R7.01.11] pour plus de détails). Les données nécessaires du champ matériau sont fournies dans l’opérateur DEFI_MATERIAU, sous les mots clés THM_LIQU et THM_GAZ.

Le tableau ci-dessous résume les mots clés obligatoires pour les sous commandes explicitées dans les sections suivantes en fonction de la loi de couplage choisie.

Légendes :

O

Mot clé Obligatoire

T

Mot clé obligatoire en Thermique

Mot clé inutile pour ce type de loi de couplage

F

Utile pour la loi mais facultatif (défini par défaut)

LIQU_SATU

LIQU_GAZ

GAZ

LIQU_GAZ_ATM

LIQU_VAPE_GAZ

LIQU_AD_GAZ_VAPE

LIQU_AD_GAZ

LIQU_VAPE

THM_INIT

O

O

O

O

O

O

O

O

PRE1

F

F

F

F

F

F

F

F

PRE2

F

F

F

F

PORO

O

O

O

O

O

O

O

O

TEMP

F

F

F

F

F

F

F

F

PRES_VAPE

O

O

O

O

THM_DIFFU

O

O

O

O

O

O

O

O

R_GAZ

O

O

O

O

O

O

RHO

O

O

O

O

O

O

O

O

BIOT_COEF

O

O

O

O

O

O

O

O

PESA_X

O

O

O

O

O

O

O

O

PESA_Y

O

O

O

O

O

O

O

O

PESA_Z

O

O

O

O

O

O

O

O

Si loi HYDR_UTIL (cf. chapitre 2.5 ) :

SATU_PRES

O

O

O

O

O

O

D_SATU_PRES

O

O

O

O

O

O

PERM_LIQU

O

O

O

O

O

O

D_PERM_LIQU_SATU

O

O

O

O

O

O

PERM_GAZ

O

O

O

O

O

D_PERM_SATU_GAZ

O

O

O

O

O

D_PERM_PRES_GAZ

O

O

O

O

O

Si loi HYDR_VGM ou HYDR_VGC (cf. section 2.5 ) :

VG_N

O

O

O

O

VG_PR

O

O

O

O

VG_SR

O

O

O

O

VG_SMAX

O

O

O

O

VG_SATUR

O

O

O

O

Si loi HYDR_ TABBAL (cf. section 2.5 ) :

SATU_PRES

O

O

O

O

O

O

D_SATU_PRES

O

O

O

O

O

O

PERM_LIQU

O

O

O

O

O

O

D_PERM_LIQU_SATU

O

O

O

O

O

O

PERM_GAZ

O

O

O

O

O

D_PERM_SATU_GAZ

O

O

O

O

O

D_PERM_PRES_GAZ

O

O

O

O

O

EPAI

O

O

O

O

A0

O

O

O

O

SHUTTLE

O

O

O

O

S_BJH

O

O

O

O

W_BJH

O

O

O

O

FICKV_T

O

O

FICKV_PV

F

F

FICKV_PG

F

F

FICKV_S

F

F

D_FV_T

F

F

D_FV_PG

F

F

FICKA_T

O

O

FICKA_PA

F

F

FICKA_PL

F

F

FICKA_S

F

F

D_FA_T

F

F

CP

T

T

T

T

T

T

T

T

PERM_IN/PERM_END/PERM_L,PERM_N, PERM_T

O

O

O

O

O

O

O

O

LAMB_T/LAMB_TL,LAMB_TN, LAMB_TT

T

T

T

T

T

T

T

T

LAMB_S

LAMB_PHI

LAMB_CT/LAMB_C_L,LAMB_C_N,LAMB_C_T

D_LB_T/D_LB_TT,D_LB_TL,D_LB_T,

D_LB_S

D_LB_PHI

THM_LIQU

O

O

O

O

O

O

O

RHO

O

O

O

O

O

O

O

UN_SUR_K

O

O

O

O

O

O

O

VISC

O

O

O

O

O

O

O

D_VISC_TEMP

O

O

O

O

O

O

O

ALPHA

T

T

T

T

T

T

T

CP

T

T

T

T

T

T

T

THM_GAZ

O

O

O

O

O

O

MASS_MOL

O

O

O

O

O

O

VISC

O

O

O

O

O

O

D_VISC_TEMP

O

O

O

O

O

O

CP

T

T

T

T

T

T

THM_VAPE_GAZ

O

O

O

MASS_MOL

O

O

O

CP

O

O

O

VISC

O

O

O

D_VISC_TEMP

O

O

O

THM_AIR_DISS

O

O

CP

O

O

COEF_HENRY

O

O

Le tableau ci-dessous résume la compatibilité entre les lois de comportement et les modélisations chosies :

MODELISATION

Loi de comportement hydraulique autorisée

*_HM(_,S,D,SI)

LIQU_SATU, GAZ, LIQU_GAZ_ATM

*_HHM(_,S,D)

LIQU_GAZ, LIQU_VAPE_GAZ

*_HH2M(S,D,_SI)

LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE

*_THH(S,D)

LIQU_GAZ, LIQU_VAPE_GAZ

*_THH2(S,D)

LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE

*_THM(_,S,D)

LIQU_SATU, GAZ, LIQU_GAZ_ATM

*_THV(_,S,D)

LIQU_VAPE

*_THHM(S,D)

LIQU_GAZ, LIQU_VAPE_GAZ

*_THH2M(S,D)

LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE

*_HH(S,D)

LIQU_GAZ, LIQU_VAPE_GAZ

*_HH2(S,D)

LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE

Mot clé facteur THM_INIT#

Pour tous les comportements Thermo-Hydro-Mécaniques, il permet de décrire un état de référence de la structure (cf. [R7.01.10] et [R7.01.11]).

Sous ce mot clé, il y a trois grandeurs correspondant à des degrés de liberté. Pour bien comprendre ces données, il faut alors distinguer les inconnues aux nœuds, que nous appelons \({\left\lbrace u\right\rbrace }^{\text{ddl}}\) et les valeurs définies sous le mot clé THM_INIT que nous appelons \({p}^{\text{ref}}\) et \({T}^{\text{ref}}\)

\({\left\lbrace u\right\rbrace }^{\text{ddl}}=\left\lbrace \begin{array}{}{u}_{x}\\ {u}_{y}\\ {u}_{z}\\ {\text{PRE1}}^{\text{ddl}}\\ {\text{PRE2}}^{\text{ddl}}\\ {T}^{\text{ddl}}\end{array}\right\rbrace\)

La signification des inconnues PRE1 et PRE2 varie suivant les modèles. En notant \({p}_{w}\) la pression d’eau, \({p}_{\text{ad}}\) la pression d’air dissous, \({p}_{\text{lq}}\) la pression de liquide \({p}_{\text{lq}}={p}_{w}+{p}_{\text{ad}}\) , \({p}_{\text{as}}\) , \({p}_{\text{vp}}\) la pression de vapeur, \({p}_{\text{as}}\) la pression d’air sec et \({p}_{g}={p}_{\text{as}}+{p}_{\text{vp}}\) la pression totale de gaz et \({p}_{c}={p}_{g}-{p}_{\text{lq}}\) la pression capillaire (aussi appelée succion), on a les significations suivantes des inconnues PRE1 et PRE2

Tableau 7 contenu de PRE1 et PRE2#

Comportement KIT

LIQU_SATU

LIQU_GAZ_ATM

GAZ

LIQU_VAPE_GAZ

\(\mathit{PRE1}\)

\({p}_{\text{lq}}\)

\(-{p}_{\text{lq}}\)

\({p}_{g}\)

\({p}_{c}={p}_{g}-{p}_{\text{lq}}\)

\(\mathit{PRE2}\)

\({p}_{g}\)
Tableau 8 contenu de PRE1 et PRE2#

Comportement KIT

LIQU_GAZ

LIQU_VAPE

LIQU_AD_GAZ_VAPE ou LIQU_AD_GAZ

\(\mathit{PRE1}\)

\({p}_{c}={p}_{g}-{p}_{\text{lq}}\)

\({p}_{\text{lq}}\)

\({p}_{c}={p}_{g}-{p}_{\text{lq}}\)

\(\mathit{PRE2}\)

\({p}_{g}\)

\({p}_{g}\)

On définit alors les pressions et la température «totales» par:

\(p={p}^{\text{ddl}}+{p}^{\text{ref}};T={T}^{\text{ddl}}+{T}^{\text{ref}}\)

Toutes les valeurs en entrée ou sortie (conditions aux limites ou résultat de IMPR_RESU ) sont les inconnues nodales \({p}^{\text{ddl}}\text{et}{T}^{\text{ddl}}\) .

Par contre ce sont les pressions et la température totale qui sont utilisées dans les lois de comportement \(\frac{p}{\rho}=\frac{R}{M}T\) pour les gaz parfaits, \(\frac{d{\rho}_{l}}{{\rho}_{l}}=\frac{{\mathrm{dp}}_{l}}{{K}_{l}}-3{\alpha}_{l}\mathrm{dT}\) pour le liquide et dans la relation saturation/pression capillaire.

On peut donc choisir d’initialiser les les valeurs nodales dans THM_INIT (par défaut ces valeurs sont nulles) ou par le mot clé ETAT_INIT de la commande STAT_NON_LINE ou par une combinaison des deux (cela peut permettre une facile conversion de Kelvin en degré Celsius par exemple).

TEMP est la température de référence \({T}_{\text{ref}}\) . Par défaut elle est prise égale à zéro. Attention la valeur de la température initiale \(T={T}^{\text{ddl}}+{T}^{\text{ref}}\) doit être strictement supérieure à zéro.

PRE1 est prise égale à zéro par défaut. Comme vu dans le tableau:

  • Pour les comportements LIQU_SATU, et LIQU_VAPE c’est la pression de liquide de référence.

  • Pour le comportement GAZ c’est la pression de gaz de référence. Dans ce cas la pression initiale de gaz :math:`p={p}^{text{ddl}}+{p}^{text{ref}}`doit être non nulle.

  • Pour le comportement LIQU_GAZ_ATM c’est la pression de liquide de référence changée de signe.

  • Pour les comportements LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE et LIQU_GAZ c’est la pression capillaire de référence.

PRE2 est prise égale à zéro par défaut. Pour les comportements LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE et LIQU_GAZ, c’est la pression de gaz de référence. La pression initiale de gaz \(p={p}^{\text{ddl}}+{p}^{\text{ref}}\) doit être non nulle.

PORO est la porosité initiale.

PRES_VAPE est la pression de vapeur de référence pour les comportements LIQU_VAPE_GAZ, LIQU_AD_GAZ_VAPE et LIQU_VAPE.

Remarques

La pression initiale de vapeur doit être prise en cohérence avec les autres données. Bien souvent, on part de la connaissance d’un état initial d’hygrométrie. Le degré hygrométrique est le rapport entre la pression de vapeur et la pression de vapeur saturante à la température considérée. On utilise alors la loi de Kelvin qui donne la pression du liquide en fonction de la pression de vapeur, de la température et de la pression de vapeur saturante:

\[\frac{{p}_{w}{-}{p}_{w}^{\text{ref}}}{{\rho}_{w}} = \frac{R}{{M}_{\mathit{vp}}^{\mathit{ol}}}T\ln(\frac{{p}_{\mathit{vp}}}{{p}_{\mathit{vp}}^{\mathit{sat}}(T)})\]

Cette relation n’est valable que pour des évolutions isothermes. On souligne que \({p}_{w}^{\text{ref}}\) correspond à un état d’équilibre auquel correspond \({p}_{\mathit{pv}}^{\mathit{sat}}\), cet état d’équilibre correspond en fait à \({p}_{w}^{0}={p}_{\mathit{gz}}^{0}=1\mathit{atm}\). Pour des évolutions avec variation de température, connaissant une loi donnant la pression de vapeur saturante à la température \({T}_{0}\), par exemple

\[{p}_{\mathit{pv}}^{\mathit{sat}}({T}_{0}) = {10}^{(2.7858+\frac{{T}_{0}{-}273.5}{31.559+0.1354({T}_{0}{-}273.5)})}\]

et un degré d’hygrométrie \(\mathit{HR}\), on en déduit la pression de vapeur grâce à \({p}_{\mathit{pv}}({T}_{0})=\mathit{HR}{p}_{\mathit{pv}}^{\mathit{sat}}({T}_{0})\). En outre, il ne faut jamais prendre une valeur de PRES_VAPE égale à zéro.

Mot clé facteur THM_LIQU#

Ce mot clé concerne tous les comportements THM faisant intervenir un liquide (cf. [R7.01.11]).

RHO est la masse volumique initiale du liquide.

UN_SUR_K est l’inverse de la compressibilité du liquide: \({K}_{l}\) .

VISC est la viscosité du liquide. Fonction de la température.

D_VISC_TEMP est la dérivée de la viscosité du liquide par rapport à la température, fonction de la température. L’utilisateur doit assurer la cohérence avec la fonction associée au mot-clef VISC.

ALPHA est le coefficient de dilatation (linéique) du liquide \({\alpha}_{l}\)

Si \({p}_{l}\) désigne la pression du liquide, \({\rho}_{l}\) sa masse volumique et \(T\) la température, le comportement du liquide est

\[\frac{{\mathrm{d\rho }}_{l}}{{\rho}_{l}}=\frac{{\text{dp}}_{l}}{{K}_{l}}-3{\alpha}_{l}\text{dT}\]

CP est la chaleur massique à pression constante du liquide.

Mot clé facteur THM_GAZ#

Ce mot clé facteur concerne tous les comportements THM faisant intervenir un gaz (cf. [R7.01.11]). Pour les comportements faisant intervenir à la fois un liquide et un gaz, et quand on prend en compte l’évaporation du liquide, les coefficients renseignés ici concernent le gaz sec. Les propriétés de la vapeur seront renseignées sous le mot clé THM_VAPE_GAZ.

MASS_MOL est la masse molaire du gaz sec. \({M}_{\text{gs}}\)

Si \({p}_{\text{gs}}\) désigne la pression du gaz sec, \({\rho}_{\mathrm{gs}}\) sa masse volumique, \(R\) la constante des gaz parfaits et \(T\) la température, le comportement du gaz sec est

\[\frac{{p}_{\text{gs}}}{{\rho}_{\text{gs}}}=\frac{\text{RT}}{{M}_{\text{gs}}}\]

CP est la chaleur massique à pression constante du gaz sec.

VISC est la viscosité du gaz sec. Fonction de la température.

D_VISC_TEMP est la dérivée par rapport à la température de la viscosité du gaz sec. Fonction de la température, l’utilisateur doit assurer la cohérence avec la fonction associée à VISC.

Mot clé facteur THM_VAPE_GAZ#

Ce mot clé facteur concerne tous les comportements THM faisant intervenir à la fois un liquide et un gaz, et prenant en compte l’évaporation du liquide (cf. [R7.01.11]). Les coefficients renseignés ici concernent la vapeur.

MASS_MOL est la masse molaire de la vapeur. \({M}_{\text{vp}}\)

CP est la chaleur massique à pression constante de la vapeur.

VISC est la viscosité de la vapeur. Fonction de la température.

D_VISC_TEMP est la dérivée par rapport à la température de la viscosité de la vapeur. Fonction de la température. L’utilisateur doit assurer la cohérence avec la fonction associée à VISC.

Mot clé facteur THM_AIR_DISS#

Ce mot clé facteur concerne le comportement THM THM_AIR_DISS prenant en compte la dissolution de l’air dans le liquide (cf. [R7.01.11]). Les coefficients renseignés ici concernent l’air dissous.

CP est la chaleur massique à pression constante de l’air dissous.

COEF_HENRY est la constante de Henry \({K}_{H}\) , permettant de relier la concentration molaire d’air dissous \({C}_{\text{ad}}^{\text{ol}}\) (\(\mathit{moles}/{m}^{3}\) ) à la pression d’air sec: \({C}_{\text{ad}}^{\text{ol}}=\frac{{p}_{\text{as}}}{{K}_{H}}\)

Remarque

La constante de Henry que nous utilisons ici s’exprime en \(\mathit{Pa.}{m}^{3}.{\mathit{mol}}^{{-}1}\). Dans la littérature il existe différentes manières d’écrire la loi de Henry. Par exemple dans les benchmarks de l’Andra, la loi de Henry est donnée par \({\omega}_{l}^{a}=\frac{{P}_{\text{as}}}{H}\frac{{M}_{\text{as}}^{\text{ol}}}{{M}_{w}}{\rho}_{w}\) avec la concentration d’air dans l’eau que l’on peut ramener à une masse volumique telle que \({\omega}_{l}^{a}={\rho}_{\text{ad}}\). \(H\) est un coefficient qui s’exprime en \(\mathit{Pa}\). Il faudra dans ces cas écrire l’équivalence \({K}_{H}=H\frac{{M}_{w}}{{\rho}_{w}}\)

Mot clé facteur THM_DIFFU#

Ce mot-clef facteur est obligatoire pour tous les comportements THM (cf. [R7.01.11]). L’utilisateur doit s’assurer de la cohérence des fonctions et de leur dérivée.

Remarque

Pour l’orthotropie, l’utilisateur a plusieurs possibilités : se placer dans un cadre isotrope classique ou se placer dans un cas d’isotropie transverse en 3D. Dans ce cas il devra renseigner un certain nombre de données matériaux (cf. ci-dessous) dans les directions L et N du repère d’orthotropie, le repère L,T étant le repère d’isotropie (cf. R4.01.02). Dans le cas 2D, la notion d’isotropie transverse perd son sens et dans ce cas seulement, on se place dans un cadre orthotrope. Pour la partie mécanique élastique on se référera à R4.01.02 et les 3 directions seront à prendre en compte. Pour les termes conductifs (hydrauliques et thermiques), seules les directions L et T du plan seront cette fois à renseigner (correspondant au repère X,Y si le repère global coïncide avec le repère local).

Attention ! L’anisotropie en mécanique n’est valable que pour un comportement mécanique élastique (non programmé pour les autres lois).

♦    R_GAZ

Constante des gaz parfaits.

♦    RHO

Pour les comportements hydrauliques masse volumique homogénéisée initiale [R7.01.11]. C’est celle qui interviendra dans l’équation d’équilibre de quantité de mouvement et dans la partie hydraulique. Attention dans les cas avec dynamique et uniquement dans ces cas, il faudra également renseigner la masse volumique dans ELAS qui sera utilisée pour les forces d’inertie (l’utilisateur devra donc prendre garde à être cohérent entre ces deux valeurs).

◊    CP

Pour les comportements thermiques, chaleur massique à contrainte constante du solide seul (des grains).

Remarque:

Attention il s’agit ici de la chaleur massique uniquement et non pas de « \(\rho {C}_{p}\) », comme c’est fait pour d’autres commandes thermiques. La masse volumique des grains est calculée dans le code à partir de la masse volumique homogénéisée [R7.01.11].

♦    BIOT_COEF

Coefficient de Biot utilisé dans le cas isotrope.

♦    BIOT_L,BIOT_N

Dans le cas isotrope transverse (3D seulement), les coefficients de Biot pour les directions L et N du repère local d’orthotropie. Dans ce cas les deux sont obligatoires.

♦    BIOT_L,BIOT_N,BIOT_T

Dans le cas orthotrope (en 2D seulement) , les coefficients de Biot pour les directions L et N T. Dans ce cas les trois sont obligatoires.

♦    PESA_X, PESA_Y, PESA_Z

Pesanteur selon x, y et z

Remarque:

La pesanteur définie ici est celle utilisée dans l’équation de Darcy uniquement. Quand il y a des calculs mécaniques, la pesanteur est également définie dans AFFE_CHAR_MECA .Cette remarque s’applique bien sur pour les trois composantes de la pesanteur.

◊    PERM_IN                 [fonction**]

Perméabilité intrinsèque fonction de la porosité dans le cas isotrope.

La perméabilité au sens classique \(K\) , dont la dimension est celle d’une vitesse se calcule de la façon suivante:

\(K=\frac{{K}_{int}{K}_{\text{rel}}}{m}{\rho}_{l}g\) où \({K}_{int}\) est la perméabilité intrinsèque, \({K}_{\text{rel}}\) la perméabilité relative, \(m\) la viscosité, \({\rho}_{l}\) la masse volumique du liquide et \(g\) l’accélération de la pesanteur. Obligatoire dans le cas isotrope.

◊    PERMIN_L,PERMIN_N [fonction**]

Dans le cas isotrope transverse 3D, la perméabilité intrinsèque pour les directions L et N du repère local d’orthotropie. Dans ce cas les deux sont obligatoires.

◊    PERMIN_L,PERMIN_T [fonction**]

Dans le cas orthotrope en 2D , la perméabilité intrinsèque pour les directions L et T du repère local (correspondant à X,Y sans rotation d’angle). Dans ce cas les deux sont obligatoires.

◊    SATU_PRES                 [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_GAZ_ATM, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), isotherme de saturation fonction de la pression capillaire. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ). On peut en outre renseigner plusieurs isothermes selon la température en utilisant alors DEFI_NAPPE (dans le cas HYDR_UTIL).

Remarque

Pour des raisons numériques, il faut éviter que la saturation atteigne la valeur 1. Aussi il est très fortement recommandé de multiplier la fonction capillaire (comprise généralement entre 0. et 1.) par 0,999.

◊    D_SATU_PRES                 [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_GAZ_ATM, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), dérivée de la saturation par rapport à la pression. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).On peut en outre renseigner plusieurs dérivées des différentes isothermes selon la température en utilisant alors DEFI_NAPPE [U4.31.03] (dans le cas HYDR_UTIL).

◊    PERM_LIQ                 [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_GAZ_ATM, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), perméabilité relative au liquide: fonction de la saturation. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).

◊    D_PERM_LIQ_SATU         [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_GAZ_ATM, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), dérivée de la Perméabilité relative au liquide par rapport à la saturation: fonction de la saturation. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).

◊    PERM_GAZ                 [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), perméabilité relative au gaz: fonction de la saturation et de la pression de gaz. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).

◊    D_PERM_SATU_GAZ         [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), dérivée de la perméabilité au gaz par rapport a la saturation: fonction de la saturation et de la pression de gaz. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).

◊    D_PERM_PRES_GAZ         [fonction**]

Pour les comportements de matériaux non saturés (LIQU_VAPE_GAZ, LIQU_VAPE, LIQU_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE), dérivée de la perméabilité au gaz par rapport a la pression de gaz: fonction de la saturation et de la pression de gaz. Uniquement pour les lois de couplage HYDR_UTIL, HYDR_TABBAL ou HYDR_ENDO (voir section 2.5 ).

◊    VG_N        [I]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydraulique est`` HYDR_VGM`` ou HYDR_VGC (voir section 2.5 ), désigne le paramètre \(N\) de la loi de Mualem Van-Genuchten servant à définir la pression capillaire et les perméabilités relatives à l’eau et au gaz.

◊    VG_PR        [R]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydraulique est HYDR_VGM ou HYDR_VGC (voir section 2.5 ), désigne le paramètre \(\mathit{Pr}\) de la loi de Mualem Van-Genuchten servant à définir la pression capillaire et les perméabilités relatives à l’eau et au gaz.

◊    VG_SR        [R]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydrauliqueestHYDR_VGM ou HYDR_VGC (voir section 2.5 ), désigne le paramètre \(\mathit{Sr}\) de saturation résiduelle de la loi de Mualem Van-Genuchten servant à définir la pression capillaire et les perméabilités relatives à l’eau et au gaz.

◊    VG_PENTR        [R]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydrauliqueestHYDR_VGM ou HYDR_VGC (voir section 2.5 ), désigne le paramètre \({P}_{e}\) de pression d’entrée de la loi de Mualem Van-Genuchten. Le milieu ne se désature que lorsque la pression capillaire est supérieure à cette valeur (0. par défaut).

◊    VG_SMAX        [R]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydraulique est HYDR_VGM ou HYDR_VGC (voir section 2.5 ), désigne la saturation maximum pour laquelle on applique la loi de Mualem Van-Genuchten. Au delà de cette saturation les courbes de Mualem-Van Genuchten sont interpolées. Cette valeur doit être très proche de 1.

◊    VG_SATUR        [R]

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ) et dans le cas où la loi hydraulique est HYDR_VGM ou HYDR_VGC (voir section 2.5 ). Au delà de la saturation définie par VG_SMAX, la saturation est multipliée par ce facteur correctif. Cette valeur doit être très proche de 1 (voir doc. [R7.01.11]).

◊    EPAI                 [fonction**]
2.5

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ, LIQU_GAZ_ATM) et dans le cas où la loi hydrauliqueest HYDR_TABBAL(voir section ). Cette fonctiondésigne l’épaisseur de la couche adsorbée \(t(\mathit{HR})`qui peut être déterminé par une loi empirique de Badmann [:ref:`2 <2>\)] (Voir Annexe).

◊ A0 [R]

2.5

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ , LIQU_AD_GAZ , LIQU_AD_GAZ_VAPE , LIQU_GAZ , LIQU_GAZ_ATM ) et dans le cas où la loi hydrauliqueest HYDR_TABBAL (voir section ). Ce paramètre désigne la surface totale des pores (ou surface spécifique du matériaux)par unité de volume. Il est déterminé par un calcul BET [4] ou BJH [3] (Voir Annexe)

◊ SHUTTLE [R]

2.5

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ, LIQU_GAZ_ATM) et dans le cas où la loi hydrauliqueest HYDR_TABBAL (voir section ). Ce paramètre désigne un terme de la loi de Shuttleworth \({\left(\frac{\partial \gamma }{\partial {\epsilon}_{s}}\right)|}_{\mu}\), il peut être recalé sur des données expérimentales de retrait de séchage (Voir Annexe).

◊ S_BJH [fonction**]

2.5

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ, LIQU_GAZ_ATM) et dans le cas où la loi hydrauliqueest HYDR_TABBAL (voir section ). Cette fonctiondésigne la fraction volumique d’eau liquide \({S}_{\mathit{BJH}}`fonction de HR, déterminée par la méthode BJH [:ref:`3 <3>\)] (Voir Annexe).

◊ W_BJH [fonction**]

2.5

Pour les comportements de matériaux non saturés à deux inconnues (LIQU_VAPE_GAZ, LIQU_AD_GAZ, LIQU_AD_GAZ_VAPE, LIQU_GAZ, LIQU_GAZ_ATM) et dans le cas où la loi hydrauliqueest HYDR_TABBAL (voir section ). Cette fonction désigne la fraction volumique de pores non saturés \({\omega}_{\mathit{BJH}}`fonction de HR, déterminée par la méthode BJH [:ref:`3 <3>\)] (Voir Annexe).

◊    FICKV_T                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, partie multiplicative du coefficient de Fick fonction de la température pour la diffusion de la vapeur dans le mélange gazeux. Le coefficient de Fick pouvant être fonction de la saturation, la température, la pression de gaz et la pression de vapeur, on le définit comme un produit de 4 fonctions: FICKV_T, FICKV_S, FICKV_PG, FICKV_VP. Seul FICKV_T est obligatoire pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE. Voir remarque en section 2.2.8 .

◊    FICKV_S                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, partie multiplicative du coefficient de Fick fonction de la saturation pour l``a diffusion de la vapeur dans le mélange gazeux. Dans le cas où on utilise cette fonction, on recommande de prendre FICKV_S(1) = 0. Voir remarque en section 2.2.8 .

◊    FICKV_PG                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, partie multiplicative du coefficient de Fick fonction de la pression de gaz pour la diffusion de la vapeur dans le mélange gazeux. Voir remarque en section 2.2.8 .

◊    FICKV_PV                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, partie multiplicative du coefficient de Fick fonction de la pression de vapeur pour la diffusion de la vapeur dans le mélange gazeux. Voir remarque en section 2.2.8 .

◊    D_FV_T                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, dérivée du coefficient FICKV_T par rapport à la température. Voir remarque en section 2.2.8 .

◊    D_FV_PG                 [fonction**]

Pour les comportements LIQU_VAPE_GAZ et LIQU_AD_GAZ_VAPE, dérivée du coefficient FICKV_PG par rapport à la pression de gaz. Voir remarque en section 2.2.8 .

◊    FICKA_T                 [fonction**]

Pour les comportement LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, partie multiplicative du coefficient de Fick fonction de la température pour la diffusion de l’air dissous dans le mélange liquide. Le coefficient de Fick pouvant être fonction de la saturation, la température, la pression d’air dissous et la pression de liquide, on le définit comme un produit de 4 fonctions: FICKA_T, FICKA_S, FICKV_PA, FICKV_PL. Dans le cas de LIQU_AD_GAZ_VAPE, seul FICKA_T est obligatoire. Voir remarque en section 2.2.8 .

◊    FICKA_S                 [fonction**]

Pour les comportement LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, partie multiplicative du coefficient de Fick fonction de la saturation pour la diffusion de l’air dissous dans le mélange liquide. Voir remarque en section 2.2.8 .

◊    FICKA_PA                 [fonction**]

Pour les comportement LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, partie multiplicative du coefficient de Fick fonction de la pression d’air dissous pour la diffusion de l’air dissous dans le mélange liquide. Voir remarque en section 2.2.8 .

◊    FICKA_PL                 [fonction**]

Pour les comportement LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, partie multiplicative du coefficient de Fick fonction de la pression de liquide pour la diffusion de l’air dissous dans le mélange liquide. Voir remarque en section 2.2.8 .

◊    D_FA_T                 [fonction**]

Pour les comportement LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, dérivée du coefficient FICKA_T par rapport à la température. Voir remarque en section 2.2.8 .

◊    LAMB_T                 [fonction**]

Partie multiplicative de la conductivité thermique du mélange dépendant de la température (Voir remarque en section 2.2.8 .). Cette opérande est obligatoire dans le cas thermique et isotrope.

◊    LAMB_TL,LAMB_TN [fonction**]

Dans le cas isotrope transverse, parties multiplicatives de la conductivité thermique du mélange dépendant de la température pour les directions L et N du repère local d’orthotropie. Ces opérandes sont obligatoires dans le cas où on a de la thermique en isotropie transverse.

◊    LAMB_TL,LAMB_TT [fonction**]

Dans le cas orthotrope 2D, parties multiplicatives de la conductivité thermique du mélange dépendant de la température pour les directions L et T. Ces opérandes sont obligatoires dans le cas où on a de la thermique en orthotropie .

◊    LAMB_S                 [fonction**]

Partie multiplicative (égale à 1 par défaut) de la conductivité thermique du mélange dépendant de la saturation (Voir remarque en section 2.2.8 ).

◊    LAMB_PHI                 [fonction**]

Partie multiplicative (égale à 1 par défaut) de la conductivité thermique du mélange dépendant de la porosité (cf. [§2.2.9]).

◊    LAMB_CT                 [fonction**]

Partie de la conductivité thermique du mélange constante et additive dans le cas isotrope (cf. [§2.2.9]). Cette constante est égale à zéro par défaut.

◊    LAMB_C_L,LAMB_C_N [fonction**]

Dans le cas isotrope transverse, parties de la conductivité thermique du mélange constantes et additives (cf. [§2.2.9]) pour les directions L et N du repère local d’orthotropie. Ces constantes sont égales à zéro par défaut.

◊    LAMB_C_L,LAMB_C_T [fonction**]

Dans le cas orthotrope 2D, parties de la conductivité thermique du mélange constantes et additives (cf. [§2.2.9]) pour les directions L et T du repère local. Ces constantes sont égales à zéro par défaut.

◊    D_LB_T                 [fonction**]

Dérivée de la partie de la conductivité thermique du mélange dépendant de la température par rapport à la température dans le cas isotrope.

◊    D_LB_TL,D_LB_TN [fonction**]

Dans le cas isotrope transverse, dérivées des parties de la conductivité thermique du mélange dépendant de la température par rapport à la température, pour les directions L et N du repère local d’orthotropie.

◊    D_LB_TL,D_LB_TT [fonction**]

Dans le cas orthotrope 2D, dérivées des parties de la conductivité thermique du mélange dépendant de la température par rapport à la température, pour les directions L et TN du repère local d’orthotropie.

◊    D_LB_S                 [fonction**]

Dérivée de la partie de la conductivité thermique du mélange dépendant de la saturation.

◊    D_LB_PHI                 [fonction**]

Dérivée de la partie de la conductivité thermique du mélange dépendant de la porosité.

◊    EMMAG                     [fonction**]

Coefficient d’emmagasinement. Ce coefficient n’est pris en compte que dans les cas des modélisations sans mécanique.

Remarque:

Attention il est important de rappeler à l’utilisateur que les paramètres BIOT_COEF et BIOT_L, BIOT_N sont incompatibles pour une même modélisation. L’utilisateur doit renseigner le paramètre BIOT_COEF s’il choisit de mener une étude dans le cas isotrope, les paramètres BIOT_L, BIOT_N s’il choisit de mener son étude en isotrope transverse, ou BIOT_L, BIOT_N et BIOT_T dans le cas orthotrope 2D. La même règle est applicable pour les paramètres PERM_IN, LAMB_T, D_LB_T et LAMB_CT. Pour ces termes conductifs, dans le cas 2D orthotrope, seules les composantes L et T sont requises.

Récapitulatif des fonctions de couplages et leur dépendance#

Les tableaux ci-dessous rappellent les différentes fonctions et leurs éventuelles dépendances et obligations.

Mot clé facteur THM_LIQU

♦

RHO

\({\rho}_{\text{lq}}^{0}\)

♦

UN_SUR_K

\(\frac{1}{{K}_{\text{lq}}}\)

◊

ALPHA

\({\alpha}_{\mathrm{lq}}\)

◊

CP

\({C}_{\text{lq}}^{p}\)

♦

VISC

\({\mu}_{\text{lq}}(T)\)

♦

D_VISC_TEMP

\(\frac{\partial {\mu}_{\text{lq}}(T)}{\partial T}\)

Mot clé facteur THM_GAZ

♦

MASS_MOL

\({M}_{\text{as}}^{\text{ol}}\)

◊

CP

\({C}_{\text{as}}^{p}\)

♦

VISC

\({\mu}_{\text{as}}(T)\)

♦

D_VISC_TEMP

\(\frac{\partial {\mu}_{\text{as}}(T)}{\partial T}\)

Mot clé facteur THM_VAPE_GAZ

♦

MASS_MOL

\({M}_{\text{VP}}^{\text{ol}}\)

♦

CP

\({C}_{\text{vp}}^{p}\)

♦

VISC

\({\mu}_{\text{vp}}(T)\)

♦

D_VISC_TEMP

\(\frac{\partial {\mu}_{\text{vp}}(T)}{\partial T}\)

Mot clé facteur THM_AIR_DISS

♦

CP

\({C}_{\text{ad}}^{p}\)

♦

COEF_HENRY

\({K}_{H}\)

Mot clé facteur THM_INIT

◊

TEMP

\({}^{\text{init}}\text{}T\)

◊

PRE1

\({}^{\text{init}}\text{}{P}_{1}\)

◊

PRE2

\({}^{\text{init}}\text{}{P}_{2}\)

♦

PORO

\({\phi }^{0}\)

◊

PRES_VAPE

\({p}_{\text{vp}}^{0}\)
Mot clé facteur THM_DIFFU#

◊

R_GAZ

\(R\)

♦

RHO

\({r}_{0}\)

◊

CP

\({C}_{\sigma}^{s}\)

◊

BIOT_COEF

\(b\)

◊

BIOT_L

\({b}_{L}\)

◊

BIOT_N

\({b}_{N}\)

◊

BIOT_T

\({b}_{T}\)

◊

SATU_PRES

\({S}_{\text{lq}}({p}_{c})\)

◊

D_SATU_PRES

\(\frac{\partial {S}_{\text{lq}}({p}_{c})}{\partial {p}_{c}}\)

♦

PESA_X

\({F}_{x}^{m}\)

♦

PESA_Y

\({F}_{y}^{m}\)

♦

PESA_Z

\({F}_{z}^{m}\)

◊

PERM_IN

\({K}^{int}(\phi )\)

◊

PERMIN_L

\({K}^{int}(\phi )\) selon L

◊

PERMIN_N

\({K}^{int}(\phi )\) selon N

◊

PERMIN_T

\({K}^{int}(\phi )\) selon T (2D)

◊

PERM_LIQU

\({k}_{\text{lq}}^{\text{rel}}({S}_{\text{lq}})\)

◊

D_PERM_LIQU_SATU

\(\frac{\partial {k}_{\text{lq}}^{\text{rel}}({S}_{\text{lq}})}{\partial {S}_{\text{lq}}}\)

◊

PERM_GAZ

\({k}_{\text{gz}}^{\text{rel}}({S}_{\text{lq}},{p}_{\text{gz}})\)

◊

D_PERM_SATU_GAZ

\(\frac{\partial {k}_{\text{gz}}^{\text{rel}}({S}_{\text{lq}},{p}_{\text{gz}})}{\partial {S}_{\text{lq}}}\)

◊

D_PERM_PRES_GAZ

\(\frac{\partial {k}_{\text{gz}}^{\text{rel}}({S}_{\text{lq}},{p}_{\text{gz}})}{\partial {p}_{\text{gz}}}\)

◊

FICKV_T

\({f}_{\text{vp}}^{T}(T)\)

◊

FICKV_S

\({f}_{\text{vp}}^{S}(S)\)

◊

FICKV_PG

\({f}_{\text{vp}}^{\text{gz}}({P}_{g})\)

◊

FICKV_PV

\({f}_{\text{vp}}^{\text{vp}}({P}_{\text{vp}})\)

◊

D_FV_T

\(\frac{\partial {f}_{\text{vp}}^{T}(T)}{\partial T}\)

◊

D_FV_PG

\(\frac{\partial {f}_{\text{vp}}^{\text{gz}}({P}_{\text{gz}})}{\partial {P}_{\text{gz}}}\)

◊

FICKA_T

\({f}_{\text{ad}}^{T}(T)\)

◊

FICKA_S

\({f}_{\text{ad}}^{S}(S)\)

◊

FICKA_PA

\({f}_{\text{ad}}^{\text{ad}}({P}_{\text{ad}})\)

◊

FICKA_PL

\({f}_{\text{ad}}^{\text{lq}}({P}_{\text{lq}})\)

◊

D_FA_T

\(\frac{\partial {f}_{\text{ad}}^{T}(T)}{\partial T}\)

◊

LAMB_T

\({\lambda}_{T}^{T}(T)\)

◊

LAMB_TL

\({\lambda}_{T}^{T}(T)\) selon L

◊

LAMB_TN

\({\lambda}_{T}^{T}(T)\) selon N

◊

LAMB_TT

\({\lambda}_{T}^{T}(T)\) selon T (2D)

◊

D_LB_T

\(\frac{\partial {\lambda}_{T}^{T}(T)}{\partial T}\)

◊

D_LB_TL

\(\frac{\partial {\lambda}_{T}^{T}(T)}{\partial T}\) selon L

◊

D_LB_TN

\(\frac{\partial {\lambda}_{T}^{T}(T)}{\partial T}\) selon N

◊

D_LB_TT

\(\frac{\partial {\lambda}_{T}^{T}(T)}{\partial T}\) selon T (2D)

◊

LAMB_PHI

\({\lambda}_{\varphi}^{T}(\varphi )\)

◊

D_LB_PHI

\(\frac{\partial {\lambda}_{\varphi}^{T}(\phi )}{\partial \phi }\)

◊

LAMB_S

\({\lambda}_{S}^{T}(S)\)

◊

D_LB_S

\(\frac{\partial {\lambda}_{S}^{T}(S)}{\partial S}\)

◊

LAMB_CT

\({\lambda}_{\text{CT}}^{T}\)

◊

LAMB_C_L

\({\lambda}_{\text{CT}}^{T}\) selon L

◊

LAMB_C_N

\({\lambda}_{\text{CT}}^{T}\) selon N

◊

LAMB_C_T

\({\lambda}_{\text{CT}}^{T}\) selon T (2D)

Remarque:

Dans le cas où il y a de la thermique:

\({\lambda}^{T}\) est fonction de la porosité, de la saturation et de la température et est donnée sous la forme du produit de trois fonctions:

\({\lambda}^{T}={\lambda}_{\varphi}^{T}(\phi ).{\lambda}_{S}^{T}({S}_{\text{lq}}).{\lambda}_{T}^{T}(T)+{\lambda}_{\text{cte}}^{T}\) avec le tenseur \({\lambda}_{T}^{T}(T)\) obligatoire et les autres fonctions par défaut prises égales à un, sauf le tenseur \({\lambda}_{\text{cte}}^{T}=0\) .

Pour le coefficient de Fick du mélange gazeux, dans le cas LIQU_VAPE_GAZet LIQU_AD_GAZ_VAPE \({F}_{\text{vp}}({P}_{\text{vp}},{P}_{\text{gz}},T,S)={f}_{\text{vp}}^{\text{vp}}({P}_{\text{vp}}).{f}_{\text{vp}}^{\text{gz}}({P}_{\text{gz}}).{f}_{\text{vp}}^{T}(T).{f}_{\text{vp}}^{S}(S)\) avec \({f}_{\text{vp}}^{T}(T)\) obligatoire, les autres fonctions étant prises par défaut égales à un, et les dérivées égales à zéro. On négligera les dérivées par rapport à la pression de vapeur et à la saturation.

Dans les cas LIQU_AD_GAZ_VAPE et LIQU_AD_GAZ, le coefficient de Fick du mélange liquide sera sous la forme : \({F}_{\text{ad}}({P}_{\text{ad}},{P}_{\text{lq}},T,S)={f}_{\text{ad}}^{\text{ad}}({P}_{\text{ad}}).{f}_{\text{ad}}^{\text{lq}}({P}_{\text{lq}}).{f}_{\text{ad}}^{T}(T).{f}_{\text{ad}}^{S}(S)\) , avec \({f}_{\text{ad}}^{T}(T)\) obligatoire, les autres fonctions étant prises par défaut égales à un, et la dérivée égale à zéro. On ne considère que la dérivée par rapport à la température (les autres sont de toutes façons prises égales à zéro).

Initialisation du calcul#

Pour définir un état initial, il faut définir un état de contraintes généralisées (aux éléments), des inconnues nodales et des variables internes.

  • Dans le mot clé THM_INIT de ** **DEFI_MATERIAU, on définit des valeurs initiales pour la porosité et la pression de vapeur (lorsque le modèle tient compte de la vapeur) mais aussi si on le souhaite de référence pour les inconnues nodales.

  • Par le mot clé ** **``DEPL`` du mot clé facteur ``ETAT_INIT`` de la commande ````STAT_NON_LINE``, on affecte le champ d’initialisation des inconnues nodales. (les valeurs physiques d’initialisations seront la somme des valeurs définies sous THM_INIT et sous ** **``ETAT_INIT``).

  • Par le mot clé ``SIGM`` du mot clé facteur ** ETAT_INIT de la commande** ** STAT_NON_LINE, on affecte le champs d’initialisation des contraintes.**

  • Par le mot clé ** VARI du mot clé facteur ** ``ETAT_INIT`` ** on affecte (éventuellement) le champs d’initialisation des variables internes.**

Afin de préciser les choses, on rappelle à quelle catégorie de variables appartient chaque grandeur physique(ces grandeurs physiques existant ou non suivant la modélisation choisie) :

Inconnues nodales

\({p}_{c},{p}_{g},{p}_{\text{lq}},T,{u}_{x},{u}_{y},{u}_{z}\)

Contraintes aux points de Gauss

\({\sigma}_{xx}^{'},{\sigma}_{yy}^{'},{\sigma}_{zz}^{'},{\sigma}_{xy}^{'},{\sigma}_{xz}^{'},{\sigma}_{yz}^{'},{\sigma}_{{p}_{xx}},{\sigma}_{{p}_{yy}},{\sigma}_{{p}_{zz}},{\sigma}_{{p}_{xy}},{\sigma}_{{p}_{xz}},{\sigma}_{{p}_{yz}},\) \(\begin{array}{}{m}_{w}{\text{,M}}_{{w}_{x}}{\text{,M}}_{{w}_{y}}{M}_{{w}_{z}},{m}_{\text{vp}}{\text{,M}}_{{\text{vp}}_{x}}{\text{,M}}_{{\text{vp}}_{y}}{M}_{{\text{vp}}_{z}},{m}_{\text{as}}{\text{,M}}_{{\text{as}}_{x}}{\text{,M}}_{{\text{as}}_{y}}{M}_{{\text{as}}_{z}},\\ {m}_{\text{ad}}{\text{,M}}_{{\text{ad}}_{x}}{\text{,M}}_{{\text{ad}}_{y}}{M}_{{\text{ad}}_{z}},{h}_{w}^{m},{h}_{\text{vp}}^{m},{h}_{\text{as}}^{m},{h}_{\text{ad}}^{m},{Q}^{'},{q}_{x},{q}_{y},{q}_{z}\end{array}\)

Variables internes

\(\phi ,{\rho}_{\text{lq}},{p}_{\text{vp}},{S}_{\text{lq}}\)

La correspondance entre nom de composante Aster et grandeur physique est explicitée en [§Annexe1].

L’initialisation des inconnues nodales ainsi que la différence entre état initial et état de référence ont été décrites et détaillées dans la section 2.2.2 . On rappelle néanmoins que \(p={p}^{\text{ddl}}+{p}^{\text{ref}}\) pour les pressions PRE1 et PRE2 et \(T={T}^{\text{ddl}}+{T}^{\text{ref}}\) pour les températures, où \({p}^{\text{ref}}\) et \({T}^{\text{ref}}\) sont définis sous le mot clé THM_INIT de la commande DEFI_MATERIAU.

Le mot clé DEPL du mot clé facteur ETAT_INIT de la commande STAT_NON_LINE définit les valeurs initiales de \({\left\lbrace u\right\rbrace }^{\text{ddl}}\) . Les valeurs initiales des masses volumiques de la vapeur et de l’air sec sont définies à partir des valeurs initiales des pressions de gaz et de vapeur. On remarque que, pour les déplacements, la décomposition \(u={u}^{\text{ddl}}+{u}^{\text{ref}}\) n’est pas faite: le mot clé THM_INIT de la commande DEFI_MATERIAU ne permet donc pas de définir des déplacements initiaux. La seule façon d’initialiser les déplacements est donc de leur donner une valeur initiale par le mot clé facteur ETAT_INIT de la commande STAT_NON_LINE.

Concernant les contraintes, les champs à renseigner sont les contraintes indiquées en annexe I suivant la modélisation choisie.

Les valeurs initiales des enthalpies, qui appartiennent aux contraintes généralisées sont définies à partir du mot clé SIGM du mot clé facteur ETAT_INIT de la commande STAT_NON_LINE. L’introduction des conditions initiales est très importante, pour les enthalpies. En pratique, on peut raisonner en considérant que l’on a trois états pour les fluides :

  • l’état courant,

  • l’état de référence: c’est celui des fluides à l’état libre. Dans cet état de référence, on peut considérer que les enthalpies sont nulles,

  • l’état initial: il doit être en équilibre thermodynamique. Pour les enthalpies de l’eau et de la vapeur, on devra prendre:

\(\begin{array}{c}{\rbrace ^{\text{init}}{h}_{w}^{m}=\frac{{p}_{w}^{\text{init}}{-}{p}_{l}^{\text{ref}}}{{\rho}_{w}}=\frac{{p}_{w}^{\text{init}}{-}{p}_{\text{atm}}}{{\rho}_{w}}\\ {\rbrace ^{\text{init}}{h}_{\text{vp}}^{m}=L({T}^{\text{init}})=\text{chaleur latente de vaporisation}\\ {\rbrace ^{\text{init}}{h}_{\text{as}}^{m}=0\\ {\rbrace ^{\text{init}}{h}_{\text{ad}}^{m}=0\end{array}\)

et avec \(L(T)=2500800-2443(T-273.15)\text{J}/\text{kg}\)

Remarque:

La pression initiale de vapeur devra être prise en cohérence avec ces choix (cf 2.2.3 ).

Concernant les contraintes mécaniques, la partition des contraintes en contraintes totale et effective s’écrit :

\(\sigma =\sigma '+{\sigma}_{p}\)

où le tenseur \(\sigma\) est la contrainte totale, c.a.d celle qui vérifie: \(\text{Div}\left(\sigma \right)+r{F}^{m}=0\)

\(\sigma '\) est la contrainte effective. Pour les lois de contraintes effectives: \(d\sigma '=f\left(d\epsilon -{\alpha}_{0}\text{dT},\alpha \right)\) , où \(\epsilon =\frac{1}{2}\left(\nabla u+{}^{T}\nabla u\right)\) et \(\alpha\) représente les variables internes.

Les composantes du tenseur \({\sigma}_{p}\) se calculent en fonction des pressions hydrauliques. L’écriture adoptée est incrémentale et, si l’on veut que les valeurs des composantes de \({\sigma}_{p}\) soient cohérentes avec les valeursinitiales de pressions définies sous les mots-clés THM_INIT et ETAT_INIT, il faut initialiser les composantes de \({\sigma}_{p}\) par le mot clé SIGM du mot clé facteur ETAT_INIT de la commande STAT_NON_LINE.

Exemple :

Les champs de déplacements initialisés dans ETAT_INIT peuvent être définis de la manière suivante :

CHAMNO=CREA_CHAMP(MAILLAGE=MAIL,

OPERATION=”AFFE”,

TYPE_CHAM=”NOEU_DEPL_R”,

AFFE=(_F(TOUT=”OUI”,

NOM_CMP=”TEMP”,

VALE=0.0,),

_F(GROUP_NO=”SURFBO”,

NOM_CMP=”PRE1”,

VALE=7.E7,),

_F(GROUP_NO=”SURFBG”,

NOM_CMP=”PRE1”,

VALE=3.E7,),

_F(GROUP_NO=”SURFBO”,

NOM_CMP=”PRE2”,

VALE=0.0,),

_F(GROUP_NO=”SURFBG”,

NOM_CMP=”PRE2”,

VALE=0.0,),),);

Et les champs de contrainte de la manière suivante :

SIGINIT=CREA_CHAMP(MAILLAGE=MAIL,

OPERATION=”AFFE”,

TYPE_CHAM=”CART_SIEF_R”,

AFFE=(_F(GROUP_MA=”BO”,

NOM_CMP=

(“SIXX”,”SIYY”,”SIZZ”,”SIXY”,”SIXZ”,

“SIYZ”,”SIPXX”,”SIPYY”,”SIPZZ”,”SIPXY”,

“SIPXZ”,”SIPYZ”,

“M11”,”FH11X”,”FH11Y”,”ENT11”,

“M12”,”FH12X”,”FH12Y”,”ENT12”,

“QPRIM”,”FHTX”,”FHTY”,”M21”,

“FH21X”,”FH21Y”,”ENT21”,

“M22”,”FH22X”, “FH22Y”,”ENT22”,),

VALE=

(0.0,0.0,0.0,0.0,0.0,

0.0,0.0,0.0,0.0,0.0,

0.0,0.0,

0.0,0.0,0.0,0.0,

0.0,0.0,0.0, 2500000.0,

0.0,0.0,0.0,0.0,0.0,0.0,0.0,

0.,0.,0.,0.),),),);

Chargements et conditions aux limites#

Toutes les conditions aux limites ou chargement sont affectées via la commande AFFE_CHAR_MECA [U4.44.01]. Les chargements sont ensuite activés par le mot clé facteur EXCIT de la commande STAT_NON_LINE.

De manière classique, trois types de conditions aux limites sont possibles:

  • Des conditions de type Dirichlet qui consistent à imposer sur une partie de frontière des valeurs fixées pour des inconnues principales appartenant à \({\left\lbrace u\right\rbrace }^{\text{ddl}}\) (et non \(u={u}^{\text{ddl}}+{u}^{\text{init}}\) ) pour cela on utilise le mot clé facteur DDL_IMPO ou FACE_IMPO de AFFE_CHAR_MECA.

  • Des conditions de type Neumann qui consistent à imposer des valeurs aux «quantités duales», soit en ne disant rien (flux nuls en hydraulique et thermique), soit en leur donnant une valeur via les mots clés FLUN, FLUN_HYDR1 et FLUN_HYDR2 du mot clé facteur FLUX_THM_REP de la commande AFFE_CHAR_MECA. Ce flux est ensuite multiplié par une fonction du temps (par défaut égale à 1.) appelée par FONC_MULT dans le sous mot clé EXCIT de la commande STAT_NON_LINE. FLUN, FLUN_HYDR1 et FLUN_HYDR2 représentent respectivement les flux thermiques, les flux d’eau et les flux de composant gazeux (cf. fin du paragraphe).

  • Des conditions de type mixtes(ou « d’échange ») qui consistent – en non saturé –à imposer une condition «extérieure» pour chaque inconnue principale d’hydraulique (\({p}_{c},{p}_{g}\) ) ainsi que les coefficients d’échange correspondant. In fine, cela consiste à imposer une valeur aux «quantités duales» fonction des pressions extérieures et des coefficients d’échange (de manière linéaire). Les pressions extérieures et les coefficients d’échanges sont renseignés respectivement par les mots clés PRE1_EXT et PRE2_EXT et COEF_11, COEF_12, COEF_21, COEF_22, en utilisant le mot clé facteur ECHANGE_THM de la commande AFFE_CHAR_MECA.

  • Les conditions mécaniques en contraintes totales \(\boldsymbol{\sigma} \, \underline{n}\) sont elles données via PRES_REP de la commande AFFE_CHAR_MECA. On se référera à la documentation de cette commande pour en connaître les possibilités.

D’un point de vue syntaxique les conditions de Dirichlet s’appliquent donc comme sur l’exemple suivant

DIRI=AFFE_CHAR_MECA(MODELE=MODELE,

DDL_IMPO=(_F(GROUP_NO=’GAUCHE’,

TEMP=0.0,),

_F(TOUT=”OUI”,

PRE2=0.0,),

_F(GROUP_NO=’GAUCHE’,

PRE1=0.0,),

_F(TOUT=”OUI”,

DX=0.0,),

_F(TOUT=”OUI”,

DY=0.0,),

_F(TOUT=”OUI”,

DZ=0.0,),

),)

Pour les conditions de Neuman, la syntaxe sera alors comme sur l’exemple suivant :

NEU1=AFFE_CHAR_MECA(MODELE=MODELE,

FLUX_THM_REP=_F(GROUP_MA=”DROIT”,

FLUN=200.,

FLUN_HYDR1=0.0,

FLUN_HYDR2=0.0),);

NEU2=AFFE_CHAR_MECA(MODELE=MODELE,

PRES_REP=_F(GROUP_MA=”DROIT”,

PRES=2.,),);

On définit ensuite la fonction multiplicativeque l’on veut appliquer, par exemple à \(\mathit{NEU1}\) :

FLUX=DEFI_FONCTION(NOM_PARA=”INST”,

VALE=

(0.0, 386.0,

315360000.0, 312.0,

9460800000.0,12.6),);

Les chargements sont ensuite activés dans STAT_NON_LINE via le mot clé EXCIT de la manière suivante:

EXCIT=(

_F(CHARGE=DIRI,),

_F(CHARGE=NEU2,),

_F(CHARGE=NEU1,

FONC_MULT=FLUX,),

),

FLUN correspond à la valeur du flux de chaleur; FLUN_HYDR1 et FLUN_HYDR2 correspondent aux valeurs des flux hydrauliques associés aux pressions \(\mathit{PRE1}\) et \(\mathit{PRE2}\). S’il n’y a pas d’ambiguïté pour la thermique ou la mécanique, en revanche les inconnues principales hydrauliques PRE1 et PRE2 changent suivant le couplage choisi. Comme on le rappelle ci-dessous :

Comportement

LIQU_SATU

LIQU_VAPE

LIQU_GAZ_ATM

GAZ

LIQU_VAPE_GAZ LIQU_GAZ LIQU_AD_GAZ_VAPE LIQU_AD_GAZ

\(\mathit{PRE1}\)

\({p}_{\text{lq}}\)

\({p}_{\text{lq}}\)

\(-{p}_{\text{lq}}\)

\({p}_{g}\)

\({p}_{c}={p}_{g}-{p}_{\text{lq}}\)

\(\mathit{PRE2}\)

\({p}_{g}\)

Les flux associés sont:

Pour PRE1, FLUN_HYDR1 : \(\left({M}_{w}+{M}_{\text{vp}}\right).n={M}_{w}^{\text{ext}}+{M}_{\text{vp}}^{\text{ext}}\)

Pour PRE2, FLUN_HYDR2 : \(\left({M}_{\text{ad}}+{M}_{\text{as}}\right).n={M}_{\text{ad}}^{\text{ext}}+{M}_{\text{as}}^{\text{ext}}\)

Pour les conditions d’échange, qui n’existe que pour les modélisations de type non saturé * HH2* , la syntaxe sera alors comme sur l’exemple suivant :

ECH=AFFE_CHAR_MECA(MODELE=MODELE,

ECHANGE_THM=_F(GROUP_MA=”DROIT”,

COEF_11=C11,COEF_12=C12,PRE1_EXT=1.E6,

COEF_21=C21,COEF_22=C22,PRE2_EXT=1.E5,),);

Les chargements sont ensuite activés dans STAT_NON_LINE via le mot clé EXCIT de la manière suivante:

EXCIT=(_F(CHARGE=ECH, TYPE_CHARGE=”SUIV”),),

Les flux associés sont alors :

Pour PRE1, \(\left({M}_{w}+{M}_{\text{vp}}\right).n=C11({p}_{1}-{p}_{{}_{1}}^{\text{ext}})-C12({p}_{2}-{p}_{{}_{2}}^{\text{ext}})\)

Pour PRE2: \(\left({M}_{\text{ad}}+{M}_{\text{as}}\right).n=C21({p}_{1}-{p}_{{}_{1}}^{\text{ext}})-C22({p}_{2}-{p}_{{}_{2}}^{\text{ext}})\)

Nous allons résumer les diverses possibilités en distinguant le cas où on impose des valeurs à PRE1 et/ou PRE2 et celui où on travaille sur des combinaisons des 2. On signale qu’on peut bien sûr avoir différents types de conditions aux limites suivant les morceaux de frontière (groupes de nœuds ou de mailles) que l’on traite. Pour un aperçu plus complet et plus détaillé de la manière dont sont traitées les conditions aux limites dans le cas non saturé, on se reportera à la note reproduite dans l’annexe 2.

  • Cas des conditions aux limites faisant intervenir les inconnues principales PRE1 et PRE2

On résume ici le cas habituel où on impose des valeurs à PRE1 et/ou PRE2.

  • Dirichlet sur PRE1 et Dirichlet sur \(\mathit{PRE2}\)

L’utilisateur impose une valeur à \(\mathit{PRE1}\) et \(\mathit{PRE2}\) ; les flux sont des résultats de calcul.

  • Dirichlet sur PRE1 et Neuman sur \(\mathit{PRE2}\)

L’utilisateur impose une valeur à PRE1et une valeur au flux associé à PRE2 en ne disant rien sur PRE2ou en donnant une valeur à FLUN_HYDR2.

  • Dirichlet sur PRE2 et Neuman sur PRE1

L’utilisateur impose une valeur à PRE2 et une valeur au flux associé à PRE1 en ne disant rien sur PRE1 ou en donnant une valeur à FLUN_HYDR1.

  • Neuman sur PRE2 et Neuman sur PRE1

Les deux flux sont imposés soit en ne disant rien sur PRE1 et/ou PRE2 (flux nuls) soit en donnant une valeur à FLUN_HYDR1 et/ou FLUN_HYDR2.

  • Cas des conditions aux limites faisant intervenir une relation linéaire entre les inconnues principales PRE1et PRE2

Il est également possible de manier des combinaisons linéaires de PRE1 et PRE2. Il faut toutefois manier cela avec précaution de manière à partir d’un problème correctement posé. La syntaxe de cet opérateur est détaillée dans la documentation de AFFE_CHAR_MECA, l’exemple ci-dessous illustre ce type de condition:

P_DDL=AFFE_CHAR_MECA(MODELE=MODELE,

LIAISON_GROUP=(_F(

GROUP_NO_1= “BORDS”,

GROUP_NO_2= “BORDS”,

DDL_1=”PRE1”,

DDL_2=”PRE2”,

COEF_MULT_1 = x,

COEF_MULT_2 = y.,

COEF_IMPO =z,),),

);

Cette commande signifie que sur la frontière définie par le groupe de nœuds \(\mathit{BORDS}\) , les pressions PRE1 et PRE2 sont reliées par la relation linéaire

\(x\mathit{PRE1}+y\mathit{PRE2}=z\)

Remarque:

Les flux imposés sont des quantités scalaires qui peuvent s’appliquer sur une ligne ou une surface interne au solide modélisé. Dans ce cas, ces conditions aux limites correspondent à une source.

Le calcul non linéaire#

La résolution est effectuée par une méthode de couplage (fiable et robuste).

Le cœur de la résolution est la commande STAT_NON_LINE. A cette commande on affecte le modèle (mot clé MODELE), les matériaux (mot clé CHAM_MATER), le/les chargements (mot clé EXCIT) et l’état initial (mot clé ETAT_INIT) que l’on a définis par toutes les commandes décrites précédemment.

Pour les informations générales concernant cette commande et sa syntaxe, on se reportera à sa documentation. On précise juste que la méthode de calculest une méthode de Newton. On ne parle donc ci-dessous que de ce qui est spécifique aux calculs THM à savoir les mots clés facteurs RELATION et RELATION_KIT``du mot clé ``COMPORTEMENT qui sont étroitement liés.

Attention :

Sous le mot clé facteur ``NEWTON``, on doit mettre une matrice de type ``TANGENTE``et non ``ELASTIQUE``.

On présente d’abord des exemples d’utilisation pour les deux stratégies de résolution.

Résolution par couplage#

Dans ce cas, la résolution est effectuée par une seule commande STAT_NON_LINE :

U0=STAT_NON_LINE(MODELE=MODELE,

CHAM_MATER=CHMAT0,

EXCIT=(

_F(CHARGE=T_IMP,),

_F(CHARGE=CALINT,

FONC_MULT=FLUX,),),

COMPORTEMENT=_F(RELATION=”KIT_THHM”,

RELATION_KIT=

(“ELAS”,”LIQU_GAZ” ,”HYDR_UTIL”),),

ETAT_INIT=_F(DEPL=CHAMNO,

SIGM=SIGINIT),

INCREMENT=_F(LIST_INST=INST1,),

NEWTON=_F(MATRICE=”TANGENTE”,REAC_ITER=1,),

CONVERGENCE=_F(RESI_GLOB_RELA=1.E-06,

ITER_GLOB_MAXI=15,

ITER_INTE_MAXI=5,),

ARCHIVAGE=_F(PAS_ARCH=1,),);

Conseils généraux d’utilisation#

RELATION est renseigné par des relations de types KIT_XXXX qui permettent de résoudre simultanément de deux à quatre équations d’équilibre. Les équations considérées dépendent du suffixe XXXX avec la règle suivante :

M désigne l’équation d’équilibre mécanique,

T désigne l’équation d’équilibre thermique,

H désigne une équation d’équilibre hydraulique,

V désigne la présence d’une phase sous forme vapeur (en plus du liquide).

Une seule lettre \(H\) signifie que le milieu poreux est saturé (une seule variable de pression \(p\) ), par exemple soit de gaz, soit de liquide, soit d’un mélange liquide/gaz (dont la pression du gaz est constante).

Deux lettres \(H\) signifient que le milieu poreux est non saturé (deux variables de pression \(p\) ), par exemple un mélange liquide/vapeur/gaz.

La présence des deux lettres \(\mathit{HV}\) signifie que le milieu poreux est saturé par un composant (en pratique de l’eau), mais que ce composant peut être sous forme liquide ou vapeur. Il n’y a alors qu’une équation de conservation de ce composant, donc un seul degré de liberté pression, mais il y a un flux liquide et un flux vapeur. Les relations possibles sont alors les suivantes :

/    'KIT_HM'

/    'KIT_THM'

/    'KIT_HHM'

/    'KIT_THH'

/    'KIT_THV'

/    'KIT_THHM'

/    'KIT_HH'

Le tableau ci-dessous résume à quel kit correspond chaque modélisation:

KIT_HM

D_PLAN_HM, D_PLAN_HMD, D_PLAN_HMS, D_PLAN_HM_SI, AXIS_HM, AXIS_HMD, AXIS_HMS, 3D_HM, 3D_HMD , 3D_HMS, 3D_HM_SI

KIT_THM

D_PLAN_THM, D_PLAN_THMD, D_PLAN_THMS, AXIS_THM, AXIS_THMD, AXIS_THMS, 3D_THM, 3D_THMD, 3D_THMS

KIT_HHM

D_PLAN_HHM, D_PLAN_HHMD, D_PLAN_HHMS, AXIS_HHM, AXIS_HHMD, AXIS_HHMS, 3D_HHM, 3D_HHMD, 3D_HHMS, D_PLAN_HH2MD, AXIS_HH2MD, 3D_HH2MD, D_PLAN_HH2MS, AXIS_HH2MS, 3D_HH2MS, D_PLAN_HH2M_SI, 3D_HH2M_SI

KIT_THH

D_PLAN_THHD, D_PLAN_THHS, AXIS_THHD, AXIS_THHS, 3D_THHD, 3D_THHS, D_PLAN_THH2D, AXIS_THH2D, 3D_THH2D, D_PLAN_THH2S, AXIS_THH2S, 3D_THH2S

KIT_THV

D_PLAN_THVD , AXIS_THVD, 3D_THVD

KIT_THHM

D_PLAN_THHMD, D_PLAN_THHMS, AXIS_THHMD, AXIS_THHMS, 3D_THHM, 3D_THHMD, 3D_THHMS, D_PLAN_THH2MD, AXIS_THH2MD, 3D_THH2MD, D_PLAN_THH2MS, AXIS_THH2MS, 3D_THH2MS

KIT_HH

D_PLAN_HHD, AXIS_HHD, 3D_HHD, D_PLAN_HH2D, AXIS_HH2D, 3D_HH2D, D_PLAN_HHS, AXIS_HHS, 3D_HHS, D_PLAN_HH2S, AXIS_HH2S, 3D_HH2S

Pour chaque phénomène modélisé (hydraulique et/ou mécanique), on doit préciser dans RELATION_KIT :

  • Le modèle de comportement mécanique du squelette s’il y a modélisation mécanique (M),

/ “ELAS’

/ “MOHR_COULOMB’

/ “CJS”

/ “HUJEUX”

/ “CAM_CLAY”

/ “BARCELONE”

/ “LAIGLE”

/ “LETK”

/ “LKR”

/ “DRUCK_PRAGER”

/ “DRUCK_PRAG_N_A”

/ “VISC_DRUCK_PRAG”

/ “GONF_ELAS”

/ “HOEK_BROWN_TOT”

/ “HOEK_BROWN_EFF”

/ “MAZARS_UNIL”

/ “ENDO_ISOT_BETON”

/ “NLH_CSRM”

  • Le comportement des liquides/gaz ,(le même que celui indiqué dans COMP_THM sous DEFI_MATERIAU, cf. [§2.2.2])

/    'LIQU_SATU’

/    'LIQU_GAZ’

/    'GAZ'

/    'LIQU_GAZ_ATM'

/    'LIQU_VAPE_GAZ'

/    'LIQU_AD_GAZ_VAPE'

/    'LIQU_AD_GAZ'

/    'LIQU_VAPE'

  • Concernant le comportement hydraulique, l’utilisateur dispose actuellement de 4 choix:

/    'HYDR_UTIL’

/    'HYDR_ENDO’

/ “HYDR_VGM”

/ “HYDR_VGC”

/ “HYDR_TABBAL”

  • Soit il choisit HYDR_UTIL : ce mot clé permet alors de renseigner la courbe de

saturation et sa dérivée en fonction de la pression capillaire ainsi que la perméabilité relative et sa dérivée en fonction de la saturation. L’utilisateur rentre ces courbes via une série de formules ensuite appelée par DEFI_MATERIAU (cf. 2.2.7 ).

  • Soit il choisit HYDR_ENDO : c’est la même chose mais permet l’utilisation des lois avec endommagement MAZARS_UNIL et ENDO_ISOT_BETON (la perméabilité est alors couplée à l’endommagement).

  • Soit il choisit HYDR_VGM ou HYDR_VGC ce qui signifie que les lois de capillarité et perméabilités obéiront au modèle de Mualem-Van-Genuchten prédéfini dans le code, autrement dit, la capillarité est donnée par une loi de Van Genuchten:

\({S}_{\mathit{we}}=\frac{S{-}{S}_{\mathit{wr}}}{1{-}{S}_{\mathit{wr}}}\) et \(m=1{-}\frac{1}{n}\)

\({S}_{\text{we}}=\frac{1}{{\left[1+{\left(\frac{{P}_{c}-{P}_{e}}{{P}_{r}}\right)}^{n}\right]}^{m}}\)

La perméabilité relative à l’eau s’exprime alors en intégrant le modèle de prédiction proposé par Mualem (1976) dans le modèle de capillarité de Van-Genuchten:

\({k}_{r}^{w}=\sqrt{{S}_{\mathit{we}}}{(1{-}{(1{-}{{S}_{\mathit{we}}}^{1/m})}^{m})}^{2}\)

La perméabilité au gaz est formulée de façon similairepar une loi de parker pour HYDR_VGM :

\({k}_{r}^{\mathit{gz}}=\sqrt{(1{-}{S}_{\mathit{we}})}{(1{-}{{S}_{\mathit{we}}}^{1/m})}^{\mathrm{2m}}\)

ou par une loi cubique pour HYDR_VGC :

\({k}_{r}^{\mathit{gz}}={(1{-}S)}^{3}\)

avec dans tous les cas les paramètres \(n\) , \(\mathit{Pr}\) , \(\mathit{Pe}\) , \(\mathit{Swr}\) , \(\mathit{Smax}\) définis dans DEFI_MATERIAU :

\(n\) correspond au paramètre VG_N

\(\mathit{Pr}\) correspond au paramètre VG_PR

\(\mathit{Pe}\) correspond au paramètreVG_PENTR (0 par défaut)

\(\mathit{Swr}\) correspond au paramètre VG_SR

A cela doivent être rajoutés deux paramètres supplémentaires correspondant à un traitement que l’on effectue sur ces courbes :

\(\mathit{Smax}\) qui correspond au paramètre VG_SMAX

\(\mathit{CSAT}\) qui correspond au paramètre VG_SATUR

Nous allons maintenant expliciter ce qu’est ce traitement et ce que sont \(\mathit{CSAT}\) et \(\mathit{Smax}\) . Pour \(S>\mathit{Smax}\) , ces courbes sont interpolées par un polynôme de degré 2, \(\mathit{C1}\) en \(\mathit{Smax}\) , de manière à éviter d’avoir à traiter des dérivées de valeurs infinies. En effet, pour \(S=1\) :

:math:`frac{partial {k}_{r}^{w}(S)}{partial S}=infty `

et dans le cas de HYDR_VGM :

:math:`frac{partial {k}_{r}^{mathit{gz}}(S)}{partial S}=infty `

Pour éviter d’avoir à traiter ce problème (qui n’a a priori pas de signification physique) on remplace ces fonctions à partir d’une saturation \(\mathit{Smax}\) par un polynôme du second ordre \(\mathit{C1}\) en ce point.

Ce qui donne pour la fonction \({k}_{r}^{w}(S)\) :

Pour \(S=\mathit{Smax}\) , on détermine le polynôme \(\mathit{PL}(S)\) tel que

\(\lbrace \begin{array}{c}\mathit{PL}({S}_{max})={k}_{r}^{w}({S}_{max})\\ \mathit{PL}'({S}_{max})=\frac{\partial {k}_{r}^{w}}{\partial S}({S}_{max})\end{array}\) et \(\mathit{PL}(1)=0\)

Pour \(S>\mathit{Smax}\) , \({k}_{r}^{w}(S)\) est remplacé par \(\mathit{PL}(S)\) .

Et, dans le cas de HYDR_VGM, pour \({k}_{r}^{\mathit{gz}}(S)\) :

Pour \(S=\mathit{Smax}\) , on détermine le polynôme \(\mathit{PG}(S)\) tel que

\(\lbrace \begin{array}{c}\mathit{PG}({S}_{max})={k}_{r}^{\mathit{gz}}({S}_{max})\\ \mathit{PG}'({S}_{max})=\frac{\partial {k}_{r}^{\mathit{gz}}}{\partial S}({S}_{max})\end{array}\) et \(\mathit{PG}(0)=1\)

Pour \(S>\mathit{Smax}\) , \({k}_{r}^{\mathit{gz}}(S)\) est remplacé par \(\mathit{PG}(S)\) comme l’exemple . Dans le cas d’une loi cubique, aucun traitement spécial n’est nécessaire.

../../../../_images/100002000000026E0000017AD1821E346362C2A5.png

Figure 2.5.2-a : Exemple de prolongement d’une perméabilité relative

Pour la succion \(S(\mathit{Pc})\) et pour \(\mathit{Pc}<\mathit{Pcmin}\) (avec \(S(\mathit{Pcmin})=\mathit{Smax}\) ) on prolonge la courbe \(S(\mathit{Pc})\) par une hyperbole telle que la courbe soit \(\mathit{C1}\) en ce point:

Pour \(S>\mathit{Smax}\) :

\(S(\mathit{Pc})=1{-}\frac{A}{B{-}{P}_{c}}\)

Avec \(A\) et \(B\) tels que la courbe soit \(\mathit{C1}\) en \(\mathit{Smax}\) .

On a donc bien une courbe décroissante qui tend vers 1 quand \(\mathit{Pc}\) tend vers :math:`{-}infty ` . Ce traitement nous permet de gérer des pressions capillaires négatives dans Code_Aster (voir l’exemple ).

../../../../_images/10000200000001DD0000012D64F26F127F8E7AED.png

Figure 2.5.2-b : Exemple de prolongement pour PC<0

\(S(\mathit{Pc})\) est ensuite multipliée par un coefficient «de sécurité» \(\mathit{CSAT}\) de manière à ce que la saturation n’atteigne jamais 1 (problème que l’on ne sait pas traiter).

Conseil : on conseille une valeur de \(\mathit{CSAT}\) le plus proche possible de \(1.\) ( \(\mathrm{0.99999}\) par exemple).

  • Soit il choisit HYDR_TABBAL qui fait appel à une formulation poromécanique valable pour des gammes d’humidités relatives au dessous de 50% (modèle issu de la thèse de Ginger El Tabbal [1]). Dans ce cas, les paramètres correspondant à la courbe de saturation et sa dérivée ainsi les perméabilitésrelativeset leur dérivée sont définis par l’utilisateur de la même manière qu’avec HYDR_UTIL. Par ailleurs, tous les paramètres spécifiques à la loi HYDR_TABBAL \(t(\mathit{HR})\) , \({A}_{0}\) , \({S}_{\mathit{BJH}}\) , \({\omega}_{\mathit{BJH}}\) , et \({\left(\frac{\partial \gamma }{\partial {\epsilon}_{s}}\right)|}_{\mu}\) doivent être définis dans DEFI_MATERIAU :

    \(t(\mathit{HR})\) correspond au paramètre EPAI

    \({A}_{0}\) correspond au paramètre A0

    \({S}_{\mathit{BJH}}\) correspond au paramètre `` S_BJH``

    \({\omega}_{\mathit{BJH}}\) correspond au paramètre W_BJH

\({\left(\frac{\partial \gamma }{\partial {\epsilon}_{s}}\right)|}_{\mu}\) correspond au paramètre SHUTTLE

Une méthode de calcul de ces paramètres est définie dans l’Annexe.

Attention :

Selon le KIT_XXXX choisi, tous les comportements ne sont pas licites (par exemple si on choisit un milieux poreux non saturé, on ne peut pas affecter un comportement de type gaz parfait). toutes les combinaisons possibles sont résumées ci-dessous

Mécanique

Couplage hydraulique

Loi hydraulique

KIT_HM KIT_THM

“ELAS’,”MOHR_COULOMB”,”CJS”, “HUJEUX”,‘LAIGLE’,‘CAM_CLAY’, ’DRUCK_PRAGER’,”LEK”, “DRUCK_PRAG_N_A’, ’VISC_DRUCK_PRAG’,”GONF_ELAS”, ‘HOEK_BROWN_TOT’, ‘HOEK_BROWN_EFF’,

‘GAZ’, ‘LIQU_SATU’, ’LIQU_GAZ_ATM’

‘HYDR_UTIL’

KIT_HM KIT_THM

‘MAZARS’, ‘ENDO_ISOT_BETON’,

‘GAZ’, ‘LIQU_SATU’, ’LIQU_GAZ_ATM’

‘HYDR_ENDO’

KIT_HHM KIT_THHM

“ELAS’,”MOHR_COULOMB”,”CJS”, “HUJEUX”‘LAIGLE’, ‘CAM_CLAY’,’DRUCK_PRAGER’, “DRUCK_PRAG_N_A’, ’VISC_DRUCK_PRAG’, “GONF_ELAS”,”LEK”, ‘HOEK_BROWN_TOT’, ‘HOEK_BROWN_EFF’,”BARCELONE”

‘LIQU_GAZ’,’LIQU_VAPE_GAZ’, ’LIQU_AD_VAPE_GAZ’ ’LIQU_AD_GAZ’

‘HYDR_UTIL’ ‘HYDR_VGM’ ‘HYDR_VGC’ ‘HYDR_TABBAL’

KIT_HHM KIT_THHM

‘MAZARS_UNIL’, ‘ENDO_ISOT_BETON’,

‘LIQU_GAZ’,’LIQU_VAPE_GAZ’, ’LIQU_AD_VAPE_GAZ’ ’LIQU_AD_GAZ’

‘HYDR_ENDO’

KIT_THH KIT_HH

‘LIQU_GAZ’,’LIQU_VAPE_GAZ’, ’LIQU_AD_VAPE_GAZ’ ’LIQU_AD_GAZ’

‘HYDR_UTIL’ ‘HYDR_VGM’ ‘HYDR_VGC’ ‘HYDR_TABBAL’

KIT_THV

’LIQU_VAPE’

‘HYDR_UTIL’

Remarque :

En cas de problème de convergence, il peut être très utile d’activer la recherche linéaire. La recherche linéaire n’améliore cependant pas systématiquement la convergence, elle est donc à manier avec précaution.

Le post-traitement#

Généralités#

Le post-traitement des données en THM ne varie pas du post-traitement Aster habituel. On rappelle juste que pour toute impression des valeurs aux nœuds qui ne sont pas les inconnues nodales, il est nécessaire de calculer ces valeurspar la commande CALC_CHAMP dont on donne un exemple ci-après.

Pour les contraintes:

U0=CALC_CHAMP(reuse =U0,

MODELE=MODELE,

CHAM_MATER=CHMAT0,

TOUT_ORDRE=”OUI”,

CONTRAINTE=(“SIGM_ELNO”),

RESULTAT=U0,);

Pour les variables internes:

U0=CALC_CHAMP(reuse =U0,

MODELE=MODELE,

CHAM_MATER=CHMAT0,

TOUT_ORDRE=”OUI”,

VARI_INTERNE=(“VARI_ELNO”),

RESULTAT=U0,);

Il faut toutefois rappeler que toutes les valeurs de déplacements en sorties correspondent à \({u}^{\text{ddl}}\) et non \(u={u}^{\text{ddl}}+{u}^{\text{ref}}\) .

Il est également important de connaître le nom des contraintes et les numéros des variables internes. Tout cela est consigné dans l’annexe I.

Ainsi l’exemple suivant permet d’imprimer la masse d’eau liquide sur le groupe de nœuds \(\mathit{HAUT}\) à tous les instants.

TAB1=POST_RELEVE_T(ACTION=_F(INTITULE=”CONT”,

GROUP_NO=(‘HAUT’),

RESULTAT=U0,

NOM_CHAM=”SIEF_ELNO”,

TOUT_ORDRE=”OUI”,

NOM_CMP=(“M11”),

OPERATION=”EXTRACTION”,),);

IMPR_TABLE(TABLE=TAB1,

FICHIER=”RESULTAT”,

FORMAT=”AGRAF”,

PAGINATION=”INST”,

NOM_PARA=(“INST”,”COOR_X”,”M11”),);

L’exemple suivant permet d’imprimer les valeurs de porosité au nœud 1 et au premier instant.

TAB2=POST_RELEVE_T(ACTION=_F(INTITULE=”DEPL”,

NOEUD=”NO1”,

RESULTAT=U0,

NOM_CHAM=”VARI_ELNO”,

NUME_ORDRE=1,

NOM_CMP=(“V2”),

OPERATION=”EXTRACTION”,),);

IMPR_TABLE(TABLE=TAB2,

FICHIER=”RESULTAT”,

FORMAT=”AGRAF”,

PAGINATION=”INST”,

NOM_PARA=(“INST”,”COOR_X”,”V2”),);

Variables internes#

Outre les généralités proposées ci-dessus, il existe désormais un moyen plus convivial d””extraire les variables internes (cf. [U4.81.01]). Pour cela des champs ont été créés dont le principe est d’extraire du champ VARI_ELGA (ou VARI_ELNO pour les champs calculés aux nœuds) la variable interne qui nous intéresse via un mot clé plus parlant que \(\mathit{V1}\) , \(\mathit{V2}\) , …

En tant que post-traitement ces champs sont calculés par CALC_CHAMP. La syntaxe à utiliser est la suivante :

  • pour un champ aux éléments cham_elem

GAMP=CALC_CHAMP(RESULTAT=U0,

OPTION='VAEX_ELNO',

NOM_VARI='GAMP');

  • pour un champ aux nœuds cham_no :

GAMP=CALC_CHAMP(reuse=GAMP,

RESULTAT=GAMP,

VARI_INTERNE='VAEX_NOEU');

Puisqu’il s’agit juste d’extraire une (et une seule!!) variable interne, les cham_elem correspondants doivent avoir été calculés au préalable :

U0=CALC_CHAMP(reuse =U0,

MODELE=MODELE,

CHAM_MATER=CHMAT0,

TOUT_ORDRE=”OUI”,

VARI_INTERNE=(“VARI_ELNO”),

RESULTAT=U0,);

La liste des différents noms symboliques des variables internes est :

“DPORO”

variation de la porosité du matériau

“DRHOLQ”

variation de la masse volumique du matériau

“DPVP”

variation de la pression de vapeur

“SATLIQ”

saturation du liquide

“EVP”

déformation plastique volumique cumulée

“IND_ETA”

Indicateur d’état mécanique

“D”

Valeur de l’endommagement

“IND_END”

Indicateur d’endommagement

“TEMP_MAX”

Température maximale

“GAMP”

Déformation déviatoire plastique cumulée

“PCR”

Pression critique

“SEUIL_HYD”

Seuil hydrique

“IND_HYD”

Indicateur d’irréversibilité hydrique

“PCOHE”

Pression de cohésion

“COMP_ROC”

Comportement de la roche

“SEUIL_ISO”

Seuil isotrope

“ANG_DEV”

Angle du seuil déviatoire

“X11”

Composantes du tenseur d’écrouissage cinématique

“X22”

Composantes du tenseur d’écrouissage cinématique

“X33”

Composantes du tenseur d’écrouissage cinématique

“X12”

Composantes du tenseur d’écrouissage cinématique

“X13”

Composantes du tenseur d’écrouissage cinématique

“X23”

Composantes du tenseur d’écrouissage cinématique

“DIST_DEV”

Distance normalisée au seuil déviatoire

“DEV_SUR_CRIT”

Rapport entre le seuil déviatoire et le seuil déviatorique critique

“DIST_ISO”

Distance normalisée au seuil isotrope

“NB_ITER”

Nombre d’itérations internes

“ARRET”

Valeur du test local d’arrêt du processus itératif

“NB_REDE”

Nombre de redécoupage local du pas de temps

“SIGNE”

Signe du produit contracté de la contrainte déviatorique par la déformation plastique déviatorique

Remarque :

Lorsque la variable à extraire ne fait pas partie des variables internes des lois concernées, une alarme est émise mais le champ est tout de même affecté à R8VIDE().

Isovaleurs#

Attention : distinction résolution linéaire et maillage quadratique

On rappelle que pour tout problème THM et afin d’éviter toute ambiguïté pour l’utilisateur, le maillage utilisé est toujours quadratique. En revanche la résolution du problème hydraulique et thermique est linéaire: les inconnues \(\mathit{PRE1}\) , \(\mathit{PRE2}\) et \(\mathit{TEMP}\) n’existent donc que sur les nœuds sommets du maillage.

A l’issue du calcul, on obtientun résultat sous cette forme :

  • Les pressions \(\mathit{PRE1}\) , \(\mathit{PRE2}\) et les températures \(\mathit{TEMP}\) sont sur les nœuds sommets des éléments.

  • Les variables internes, les contraintes sont sur les points de gauss et/ou les nœuds du maillage quadratique (suivant que l’on soit en sélectif ou en lumpé) et les déplacements mécaniques sont sur les sommets et les nœuds milieux.

Il y a donc un conflit entre les maillages quadratiques et les calculs linéaires ce qui va poser un problème par exemple pour tracer des isovaleurs.

Il existe alors une manière simple de s’affranchir de ce problème en projetant l’ensemble des résultats sur le maillage linéaire issu du maillage quadratique. Ainsi tous les résultats seront définis au même endroit (sur les sommets).

Pour \(\mathit{TEMP}\) , \(\mathit{PRE1}\) et \(\mathit{PRE2}\) , il suffit de projeter comme sur l’exemple suivant:

## projection du maillage quadratique sur le maillage linéaire :

MAILQ4 = CREA_MAILLAGE(

MAILLAGE = MAIL,

QUAD_LINE = _F(TOUT = “OUI”),

)

## affectation du modèle correspondant :

MODELQ4=AFFE_MODELE(MAILLAGE=MAILQ4,

AFFE=_F(TOUT=”OUI”,

PHENOMENE=”MECANIQUE”,

MODELISATION=”D_PLAN”,),);

## projection du résultat sur ce nouveau modèle :

RESULIN = PROJ_CHAMP(

METHODE=”ELEM”,

RESULTAT=resu,

MODELE_1=MODELQ8,

MODELE_2=MODELQ4,

NOM_CHAM =”DEPL”,

TOUT_ORDRE=”OUI”)

## post-traitement habituel mais sur le maillage linéaire :

IMPR_RESU(FORMAT=”MED”,

RESU=_F(

RESULTAT=RESULIN,

NUME_ORDRE=(1,10,20,),

NOM_CMP = (“PRE1”,”PRE2”,),

MAILLAGE=MAILQ4,),);

Il existe une méthode alternative, qui peut s’avérer plus rapide pour les grosses études industrielles. Cette méthode utilise MACR_ADAP_MAIL [U7.03.01]. On suppose que resu est composé de nbinst instants stockés dans la liste linst. Le post-traitement précédent peut également se faire de la manière suivante :

PRE = [None]*nbinst;

PREB = [None]*nbinst;

for k in range(nbinst) :

PRE[k]=CREA_CHAMP(TYPE_CHAM=”NOEU_DEPL_R”,

OPERATION=”EXTR”,

RESULTAT=resu,

NOM_CHAM=”DEPL”,

INST=linst[k],);

motcle1={}

motcle2={}

motcle1[“MAJ_CHAM”]=[]

motcle2[“AFFE”]=[]

for k in range(nbinst) :

PREB[k] =CO(“PREB_%d” % (k))

motcle1[“MAJ_CHAM”].append(

_F(TYPE_CHAM=”NOEU_DEPL_R”,CHAM_GD = PRE[k],

CHAM_MAJ = PREB[k],),);

motcle2[“AFFE”].append(

_F(CHAM_GD=PREB[k],INST=linst[k],),);

MACR_ADAP_MAIL( ADAPTATION=”MODIFICATION”,DEGRE=”OUI”,

MAILLAGE_N = MAILQ8,

MAILLAGE_NP1 = CO(“MAILQ4”),

motcle1);

RESULINB=CREA_RESU(OPERATION=”AFFE”,TYPE_RESU=”EVOL_NOLI”,NOM_CHAM=”DEPL”,**motcle2);

Quelques cas tests#

A titre informatif, on récapitule ici quelques cas tests disponibles dans la base de Code_Aster. Cette liste n’est absolument pas exhaustive mais permet d’accéder directement à un certain nombre d’exemples représentatifs.

Nom de cas-test

Type de modélisation THM

Phénomène représenté

Modélisation géométrique

Loi de couplage

Wtna109a

AXIS_THVD

Désaturation par chauffage d’un milieu sans air

axisymétrique

LIQU_VAPE

Wtna110a

AXIS_HH2MS

Modélisation d’un gonflement d’une argile

axisymétrique

LIQU_AD_GAZ_VAPE

Wtna107

AXIS_THMS

Chauffage d’un milieu saturé (THM)

axisymétrique

LIQU_SATU

wtna112

AXIS_THMS

Pressurisation thermique d’une éprouvette drainée

axisymétrique

LIQU_SATU

Wtnp117a

D_PLAN_HHS

Rééquilibrage capillaire de 2 matériaux (BO/BG)

plane

LIQU_GAZ

Wtnp118a

D_PLAN_HH2D

Rééquilibrage gravitaire d’une colonne

plane

LIQU_AD_GAZ_VAPE

Wtnp112d

D_PLAN_HMD

Resaturation d’une colonne (hypothèse de Richards)

plane

LIQU_GAZ_ATM

Wtnp101a

D_PLAN_HM

Modélisation d’un milieu saturé en gaz

plane

GAZ

Wtnl100a

D_PLAN_HMD D_PLAN_HM_SI

Problème de consolidation (Therzagi)

plane

LIQU_SATU

Wtnv123a

3D_HHM

Essai hydrostatique puis triaxial à succion fixe (modèle de Barcelone)

3D

LIQU_GAZ

Wtnp113a

D_PLAN_THHD

Resaturation d’une alvéole par le massif

plane

LIQU_GAZ

Wtnv121a

3D_HHMD

Mouillage du béton avec endommagement

3D

LIQU_GAZ

Wtnv130a

3D_THH2S

Chauffage d’un milieu poreux avec air dissous

3D

LIQU_AD_GAZ_VAPE

wtnv140

3D_HMS

Essai triaxial drainé anisotrope

3D

LIQU_SATU

Les modélisations THM par une approche second gradient#

Les modélisations de type second gradient et second gradient de dilatation sont décrites dans la documentation [R5.04.03] dans le cadre des milieux poreux. Leur utilisation est indispensable lorsque les lois de comportements mécaniques modélisent une dégradation du matériau – on parle de lois adoucissantes. En effet, dans ce cas, les solutions dépendent du maillage y compris lorsque le milieu est couplé avec les équations de l’hydraulique.

L’intégration numérique de ces modèles a suivi un protocole atypique lors de sa mise en œuvre pour faciliter son utilisation dans Code_Aster. Ainsi, la partie régularisante (second gradient) est définie comme un patch qu’il suffit de «coller» sur la structure à modéliser. On décrit ci-dessous les étapes à suivre pour produire une telle simulation.

  1. Définition du patch

  2. Choix du modèle

  3. Définition du matériau

  4. Conditions aux limites

  5. Calcul

Définition du patch#

Pour la simulation d’une étude THM par une approche de type second gradient, il faut définir un maillage quadratique pour la discrétisation de la structure, puis un second maillage pour le patch régularisant. La particularité de ce second maillage est qu’il doit nécessairement être topologiquement identique au premier maillage (mêmes nœuds et mêmes éléments géométriques. Pour des raisons de simplification numérique et pour des raisons de cohérence de modèles, on conseille d’utiliser la procédure suivante:

  1. Définition d’une discrétisation spatiale pour la structure et la modélisation THM;

  2. Duplication du maillage précédent sur l’ensemble de la structure par la commande CREA_MAILLAGE pour définir le patch

  3. Modification éventuelle du type de mailles pour une étude avec multiplicateurs de Lagrange aux centres des éléments.

Étape 1. Définition du maillage de la structure#

Cette première étape est à produire par un mailleur pour discrétiser la structure. On récupère ensuite le maillage par la commande Code_Aster :

lecture du maillage (quadratique)

MA = LIRE_MAILLAGE()

Étape 2. Duplication du maillage pour définir le patch#

Il existe dans Code_Aster une procédure qui permet de dupliquer un maillage en conservant les mêmes nœuds que le maillage original – il est indispensable que les nœuds du premier maillage soient inclus dans le patch (cf [R5.04.03]). Il est nécessaire de définir ce second maillage pour permettre de définir la modélisation second gradient au patch régularisant.

Duplication du maillage (quadratique). Seules les mailles sont dupliquées, les nœuds restent communs. Le but est d’associer à chacun de ces maillages une modélisation différente.

MAIL=CREA_MAILLAGE(MAILLAGE = MA,

CREA_GROUP_MA = ( _F( NOM = “ROCHE_REG”,

GROUP_MA = “ROCHE”)))

Étape 3. Modification (éventuelle) du maillage du patch#

Il existe différentes discrétisations éléments finis pour interpoler les modélisations second gradient (avec ou sans multiplicateurs de Lagrange, multiplicateurs de Lagrange au centre des éléments ou aux nœuds sommets). Par conséquent il faut adapter le maillage du patch à l’espace de discrétisation de la modélisation second gradient.

Par exemple pour prendre en compte une interpolation avec des multiplicateurs de Lagrange aux centres des éléments, il faut définir les nœuds aux centres des éléments car ceux-ci ne sont pas prévus, en général, dans les discrétisations spatiales proposées par la plupart des mailleurs. Ainsi dans ce cas, on peut utiliser la commande Code_Aster suivante:

Introduction de nœuds centraux aux éléments finis du nouveau maillage pour une interpolation \(\mathit{P2}-\mathit{P1}-\mathit{P0}\) . Procédure nécessaire pour prendre en compte les interpolations des multiplicateurs de Lagrange.

MAILLAGE=CREA_MAILLAGE(MAILLAGE = MAIL,

MODI_MAILLE =( _F( GROUP_MA = “ROCHE_REG”,

OPTION = “QUAD8_9”),

_F( GROUP_MA = “ROCHE_REG”,

OPTION = “TRIA6_7”)))

Cette étape n’est possible qu’en 2D lorsque l’interpolation élément fini considère des multiplicateurs de Lagrange constant par élément.

Choix du modèle#

Les «patchs régularisants» second gradient ou second gradient de dilatation sont à combiner avec n’importe quel type de modélisation THM. Ils n’ont d’intérêt, bien entendu, qu’en présence d’une composante mécanique (adoucissante). Les modélisations disponibles sont:

Modélisation

Modélisation géométrique

Phénomène pris en compte

D_PLAN_2DG

Plane

Second gradient

D_PLAN_DIL

Plane

Second gradient de dilatation

3D_DIL

3D

Second gradient de dilatation

Remarque:

A noter qu’il est fortement conseillé d’utiliser les modélisations second gradient de dilatation pour traiter les comportements des sols et des roches ou plus particulièrement tout type de matériau présentant une variation volumique lors de sa dégradation .

Les inconnues nodales des modélisations THM sont alors enrichies après ajout du «patch» second gradient. Le nombre d’inconnues varie selon le patch considéré et le type de mailles caractérisant la discrétisation spatiale du patch. Pour la modélisation second gradient de dilatation c’est donc par le choix de la topologie de la maille que l’interpolation des multiplicateurs de Lagrange sera déterminée (voir le tableau ci-dessous). On résume dans le tableau ci-dessous les différentes possibilités où l’on note en accord avec [R5.04.03], \({f}_{ij}\) les composantes du tenseur des déformations microscopiques, \(\chi\) la déformation volumique microscopique et \(\lambda\) le multiplicateur de Lagrange:

Modélisation

Type de mailles

Degrés de liberté

Position du degré de liberté

D_PLAN_2DG

TRIA7, QUAD9

\({f}_{ij}\)

Aux sommets de chaque élément

\({\lambda}_{ij}\)

Au centre de chaque élément

D_PLAN_DIL

TRIA7, QUAD9

\(\chi\)

Aux sommets de chaque élément

\(\lambda\)

Au centre de chaque élément

TRIA6, QUAD8

\(\chi\)

Aux sommets de chaque élément

3D_DIL

TETRA10, HEXA20, PENTA15

\(\chi\) , \(\lambda\)

Aux sommets de chaque élément

On procède alors de la façon suivante:

MODELE = AFFE_MODELE( MAILLAGE = MAILLAGE,

AFFE = ( _F( GROUP_MA = “ROCHE”,

PHENOMENE = “MECANIQUE”,

MODELISATION = “D_PLAN_HMS”),

_F( GROUP_MA = “ROCHE_REG”,

PHENOMENE = “MECANIQUE”,

MODELISATION = “D_PLAN_DIL”)))

Définition du matériau#

Il faut, en plus de la définition des paramètres matériaux nécessaires pour le calcul THM décrit au chapitre précédent, ajouter un paramètre qui caractérise la longueur caractéristique de la modélisation second gradient (cf [R5.04.03]). Cette donnée caractérise en fait la dimension du voisinage qui contribue à la description d’un point matériel.

Cette longueur est définie à partir des composantes matériau \(\mathit{A1}\) et \(\mathit{A3}\) de la loi de comportement d’élasticité linéaire second gradient inspirée des travaux de Mindlin dans le cas de la modélisation second gradient de dilatation. C’est dans la détermination de ces paramètres que réside la plus grande difficulté. Il n’existe pas aujourd’hui de méthode analytique pour les identifier. La pratique est de faire quelques essais préalable. Dans le cas des modélisations second gradient cette longueur est fonction des 5 composantes \(\mathit{A1}\) , \(\mathit{A2}\) , \(\mathit{A3},\) \(\mathit{A4}\) et \(\mathit{A5}\) , à renseigner sous le mot-clé ELAS_2NDG de DEFI_MATERIAU (cf [R5.04.03]).

On peut, de plus, ajouter un paramètre numérique de pénalisation. Il a été montré que le terme de pénalisation des modélisations second gradient de dilatation améliore la convergence numérique sans perturber la qualité des résultats.

Voici un exemple de définition des paramètres matériau pour la partie second gradient.

SOL2 = DEFI_MATERIAU( ELAS_2NDG = _F( A1 = 4.0E4,

A2 = 0.0,

A3 = 0.0,

A4 = 0.0,

A5 = 0.0,),

NON_LOCAL = _F(PENA_LAGR=1.0E12,),);

Puis affectation des paramètres matériaux suivant la même procédure que pour la définition des modélisations.

MATE = AFFE_MATERIAU( MAILLAGE = MAILLAGE,

AFFE = (_F( TOUT = ‘ROCHE’,

MATER = SOL_HM,),

_F( GROUP_MA = “ROCHE_REG”,

MATER = SOL2)))

Impact sur les conditions aux limites#

Il n’y a aucun impact à apporter dans le fichier de commande Code_Aster pour prendre en compte les conditions aux limites avec les modélisations second gradient par rapport aux modélisations classiques THM. En revanche, il est essentiel de noter que la signification physique des conditions aux limites de pression de contraintes normales est modifiée avec les modélisations second gradient (cf [R5.04.03]).

Résolution du problème#

Le calcul est à effectuer par la commande STAT_NON_LINE. Il suffit d’ajouter au calcul classique THM la relation de comportement d’élasticité second gradient sur le maillage correspond au patch comme sur l’exemple ci-dessous (mot-clé ELAS). Il n’y a aucune incompatibilité avec les lois de comportement. A noter, tout de même, qu’une simulation THM par une approche second gradient n’a d’intérêt qu’en présence d’une composante mécanique dans le calcul couplé.

Il est important de noter que lorsque la modélisation prend en compte une interpolation des multiplicateurs de Lagrange constants par éléments il est nécessaire d’utiliser le solveur Mumps pour la résolution.

Définition du calcul statique non linéaire avec une loi de comportement associé à chacune des modélisation: comportement de type Drucker-Prager pour le premier gradient, et élasticité linéaire pour le second gradient

U1 = STAT_NON_LINE( MODELE = MODELE,

CHAM_MATER = MATE,

EXCIT = ( _F( CHARGE = CHCI ),…),

SOLVEUR = ( _F( METHODE=”MUMPS”,)),

COMPORTEMENT = ( _F(GROUP_MA=”ROCHE”,

RELATION = “DRUCK_PRAGER”,),

_F(GROUP_MA = “ROCHE_REG”,

RELATION = “ELAS”,),),

NEWTON = _F( MATRICE = “TANGENTE”,REAC_ITER = 1 ),

INCREMENT = _F( LIST_INST = TEMPS ))

On trouvera un exemple de modélisation par second gradient, par exemple dans le cas test ssll117.

Bibliographie#

    1. El TABBAL «Drying behavior of concrete: Application to nuclear waste packages». Thèse Université Paris-Est. (2020)

    1. Badmann, N. Stockausen, M.J. Setzer. «The statistical thickness and the chemical potential of adsorbed water films». J. Colloid Interface Sci 82, 534-543. 1981

    1. Barrett, L. Joyner, P. Halenda. «The determination of pore volume and area distributions in porous substances.» . Computations from Nitrogen Isotherms. 1896 . 1951

    1. Brunauer, P.H. Emett, E. Teller. «Adsorption of gases in multimolecular layers.» J. Am. Chem. Soc 60, 309-319 (1938)

  2. Rouquerol, F. et al., 2003. Texture des matériaux pulvérulents ou poreux. Tech. l’Ingénieur P1050, 1–24.

  3. G.El Tabbal «Modelling the Drying Shrinkage of Porous Materials by Considering both Capillary and Adsorption Effects. Journal of the Mechanics and Physics of Solids, Volume 142, September 2020.

Contraintes généralisées et variables internes

Les contraintes:

Numéro

Nom de composante Aster

Contenu

Modélisations

1

SIXX

\({\sigma}_{xx}^{'}\)

Si mécanique (..M…)

2

SIYY

\({\sigma}_{yy}^{'}\)

Si mécanique (..M…)

3

SIZZ

\({\sigma}_{zz}^{'}\)

Si mécanique (..M…)

4

SIXY

\({\sigma}_{xy}^{'}\)

Si mécanique (..M…)

5

SIXZ

\({\sigma}_{xz}^{'}\)

Si mécanique (..M…)

6

SIYZ

\({\sigma}_{yz}^{'}\)

Si mécanique (..M…)

7

SIPXX

\({\sigma}_{{p}_{xx}}\)

Si mécanique (..M…)

8

SIPYY

\({\sigma}_{{p}_{yy}}\)

Si mécanique (..M…)

9

SIPZZ

\({\sigma}_{{p}_{zz}}\)

Si mécanique (..M…)

10

SIPXY

\({\sigma}_{{p}_{xy}}\)

Si mécanique (..M…)

11

SIPXZ

\({\sigma}_{{p}_{xz}}\)

Si mécanique (..M…)

12

SIPYZ

\({\sigma}_{{p}_{yz}}\)

Si mécanique (..M…)

13

M11

\({m}_{w}\)

Dans tous les cas

14

FH11X

\({M}_{{w}_{x}}\)

Dans tous les cas

15

FH11Y

\({M}_{{w}_{y}}\)

Dans tous les cas

16

FH11Z

\({M}_{{w}_{z}}\)

Dans tous les cas

17

ENT11

\({h}_{w}^{m}\)

Dans tous les cas avec thermique

18

M12

\({m}_{{}_{\text{vp}}}\)

Si 2 pressions inconnues (..HH…)

19

FH12X

\({M}_{{\text{vp}}_{x}}\)

Si 2 pressions inconnues (..HH…)

20

FH12Y

\({M}_{{\text{vp}}_{y}}\)

Si 2 pressions inconnues (..HH…)

21

FH12Z

\({M}_{{\text{vp}}_{z}}\)

Si 2 pressions inconnues (..HH…)

22

ENT12

\({h}_{\text{vp}}^{m}\)

Si 2 pressions inconnues et thermique (..THH…)

23

M21

\({m}_{{}_{\text{as}}}\)

Si 2 pressions inconnues (..HH…)

24

FH21X

\({M}_{{\text{as}}_{x}}\)

Si 2 pressions inconnues (..HH…)

25

FH21Y

\({M}_{{\text{as}}_{y}}\)

Si 2 pressions inconnues (..HH…)

26

FH21Z

\({M}_{{\text{as}}_{z}}\)

Si 2 pressions inconnues (..HH…)

27

ENT21

\({h}_{\text{as}}^{m}\)

Si 2 pressions inconnues et thermique (..THH…)

28

M22

\({m}_{{}_{\text{ad}}}\)

Si modélisation de l’air dissous (…HH2…)

29

FH22X

\({M}_{{\text{ad}}_{x}}\)

Si modélisation de l’air dissous (…HH2…)

30

FH22Y

\({M}_{{\text{ad}}_{y}}\)

Si modélisation de l’air dissous (…HH2…)

31

FH22Z

\({M}_{{\text{ad}}_{z}}\)

Si modélisation de l’air dissous (…HH2…)

32

ENT22

\({h}_{\text{ad}}^{m}\)

Si modélisation de l’air dissous et thermique (…THH2…)

33

QPRIM

\(Q'\)

Si thermique

34

FHTX

\({q}_{x}\)

Si thermique

35

FHTY

\({q}_{y}\)

Si thermique

36

FHTZ

\({q}_{z}\)

Si thermique

Dans le cas sans mécanique, et pour les lois de comportements (LIQU_VAPE_GAZ,LIQU_VAPE, LIQU_AD_GAZ_VAPE , ’LIQU_AD_GAZ’) les variables internessont :

Numéro

Nom de composante Aster

Contenu

1

\(\mathit{V1}\)

\({\rho}_{\text{lq}}-{\rho}_{{\text{lq}}_{}}^{0}\)

2

\(\mathit{V2}\)

\(\varphi -{\varphi}^{0}\)

3

\(\mathit{V3}\)

\({p}_{\text{vp}}-{p}_{\text{vp}}^{0}\)

4

\(\mathit{V4}\)

\({S}_{\text{lq}}\)

Dans le cas sans mécanique, et pour les lois de comportements (LIQU_GAZ,LIQU_GAZ_ATM,) les variables internessont :

Numéro

Nom de composante Aster

Contenu

1

\(\mathit{V1}\)

\({\rho}_{\text{lq}}-{\rho}_{{}_{\text{lq}}}^{0}\)

2

\(\mathit{V2}\)

\(\varphi -{\varphi}^{0}\)

3

\(\mathit{V3}\)

\({S}_{\text{lq}}\)

Dans le cas sans mécanique, et pour les lois de comportements (LIQU_SATU,) les variables internessont :

Numéro

Nom de composante Aster

Contenu

1

\(\mathit{V1}\)

\({\rho}_{\text{lq}}-{\rho}_{{}_{\text{lq}}}^{0}\)

2

\(\mathit{V2}\)

\(\varphi -{\varphi}^{0}\)

Dans le cas avec mécanique les premiers numéros seront ceux correspondant à la mécanique (\(\mathit{V1}\) dans le cas élastique, \(\mathit{V1}\) et suivantes pour des modèles plastiques). Le numéro des variables internes ci-dessus devra alors être incrémenté d’autant.

Éléments supplémentaires sur les conditions aux limites en THM

Dans ce qui suit on ne prend pas en compte l’air dissous (l’indice \(\mathit{lq}\) correspond alors à celui de l’eau \(w\) ) et on s’attache au cas non saturé.

Nous rappelons ici le choix des inconnues de pression.

Comportement

LIQU_GAZ et LIQU_VAPE_GAZ

\(\mathit{PRE1}\)

Pression capillaire: \({p}_{c}={p}_{\text{gz}}-{p}_{\text{lq}}\)

\(\mathit{PRE2}\)

Pression de gaz \({p}_{\text{gz}}={p}_{\text{vp}}+{p}_{\text{as}}\)

A2.1 Formulation variationnelle des équations de conservation

On se réfère ici à [R7.01.11]. Ces équations sont

\(\dot{{m}_{\text{lq}}}+{\dot{m}}_{\text{vp}}+\text{Div}({M}_{\text{lq}}+{M}_{\text{vp}})=0\) éq A5.1-1

\(\dot{{m}_{\text{as}}}+\text{Div}({M}_{\text{as}})=0\) éq A5.1-2

La formulation variationnelle déduite est donnée par

\(\begin{array}{}-{\int}_{\Omega}(\dot{{m}_{\text{lq}}}+{\dot{m}}_{\text{vp}}){\pi}_{1}d\Omega +{\int}_{\Omega}({M}_{\text{lq}}+{M}_{\text{vp}}).\nabla {\pi}_{1}d\Omega =\\ {\int}_{\partial \Omega }({M}_{{\text{lq}}_{\text{ext}}}+{M}_{{\text{vp}}_{\text{ext}}}).{\pi}_{1}d\Gamma \forall {\pi}_{1}\in {P}_{{1}_{\text{ad}}}\end{array}\) éq A5.1-3

\(\begin{array}{}-{\int}_{\Omega}\dot{{m}_{\text{as}}}{\pi}_{2}d\Omega +{\int}_{\Omega}{M}_{\text{as}}.\nabla {\pi}_{2}d\Omega =\\ {\int}_{\partial \Omega }{M}_{{\text{as}}_{\text{ext}}}.{\pi}_{2}d\Gamma \forall {\pi}_{2}\in {P}_{{2}_{\text{ad}}}\end{array}\) éq A5.1-4

Les pressions capillaires et de gaz sont liées aux pression d’eau, de vapeur et d’air sec par les relations:

\({p}_{c}={p}_{\text{gz}}-{p}_{\text{lq}}\) éq A5.1-5

\({p}_{\text{gz}}={p}_{\text{vp}}+{p}_{\text{as}}\) éq A5.1-6

La pression de vapeur n’est pas une variable indépendante. Elle est reliée à la pression de liquide \({p}_{\text{lq}}\) par les relations

\(\frac{{\text{dp}}_{\text{vp}}}{{\rho}_{\text{vp}}}=\frac{{\text{dp}}_{\text{lq}}}{{\rho}_{\text{lq}}}+({h}_{\text{vp}}^{m}-{h}_{\text{lq}}^{m})\frac{\text{dT}}{T}\) éq A5.1-7

\({\text{dh}}_{\text{lq}}^{m}={C}_{\text{lq}}^{p}\text{dT}+(1-3{\alpha}_{\text{lq}}T)\frac{{\text{dp}}_{\text{lq}}}{{\rho}_{\text{lq}}}\) éq A5.1-8

\({\text{dh}}_{\text{vp}}^{m}={C}_{\text{vp}}^{p}\text{dT}\) éq A5.1-9

Ces relations montrent que la pression de vapeur est déterminée complètement pas la connaissance de \({p}_{\text{lq}}\) (et de son évolution). Souvent, ces relations servent à établir la loi de Kelvin, \(\frac{{p}_{\text{lq}}-{p}_{\text{lq}}^{0}}{{\rho}_{\text{lq}}}=\frac{R}{{M}_{\text{vp}}^{\text{ol}}}T\ln(\frac{{p}_{\text{vp}}}{{p}_{\text{vp}}^{\text{sat}}(T)})\) , mais cette loi n’est pas utilisée directement dans Aster.

Les documents de référence Aster ne disent rien sur ce que sont les variables \({\pi}_{1}\) et \({\pi}_{2}\) . Mais deux éléments peuvent nous mettre sur la piste:

  • D’une part, \({\pi}_{1}\in {P}_{{1}_{\text{ad}}}\) et \({\pi}_{2}\in {P}_{{2}_{\text{ad}}}\) alors que \({P}_{{1}_{\text{ad}}}\) et \({P}_{{2}_{\text{ad}}}\) sont les espaces d’appartenance de \(\mathit{PRE1}\) et \(\mathit{PRE2}\) (incluant donc leurs conditions aux limites).

  • D’autre part, au chapitre 7. de [R7.01.10], on voit que la déformation virtuelle \({E}_{g}^{\text{*}\mathrm{el}}=(v,\varepsilon (v),{\pi}_{1},\nabla {\pi}_{1},{\pi}_{2},\nabla {\pi}_{2},\tau ,\nabla \tau )\) est liée au vecteur de déplacement nodaux virtuel \({U}^{\text{*}\mathrm{el}}=(v,{\pi}_{1},{\pi}_{2},\tau )\) par le même opérateur \({Q}_{g}^{\text{el}}\) que celui qui relie entre eux la déformation \({E}_{g}^{\text{el}}=(u,\varepsilon (u),{p}_{1},\nabla {p}_{1},{p}_{2},\nabla {p}_{2},T,\nabla T)\) et le déplacement nodal \({U}^{\text{el}}=(u,{p}_{1},{p}_{2},T)\) :

  • \({E}_{g}^{\text{*}\mathrm{el}}={Q}_{g}^{\text{el}}{U}^{\text{*}\mathrm{el}}\)

  • \({E}_{g}^{\text{el}}={Q}_{g}^{\text{el}}{U}^{\text{el}}\)

Il est alors clair que \({\pi}_{1}\) et \({\pi}_{2}\) sont des variations virtuelles de \({p}_{1}\) et \({p}_{2}\)

D’où le tableau:

\({p}_{1}={p}_{c}\Rightarrow {\pi}_{1}=\delta {p}_{c}={p}_{c}^{\text{*}}\)

\({p}_{1}={p}_{\text{lq}}\Rightarrow {\pi}_{1}=\delta {p}_{\text{lq}}={p}_{\mathrm{lq}}^{\text{*}}\)

\({p}_{2}={p}_{\text{gz}}\Rightarrow {\pi}_{2}=\delta {p}_{\text{gz}}={p}_{\mathrm{gz}}^{\text{*}}\)

A2.2 Cas de conditions aux limites faisant intervenir les inconnues principales

Ce que nous disons dans ce paragraphe et les suivants est relatif à une partie de la frontière \(\partial {O}_{d}\) sur laquelle des conditions sont prescrites: rien n’empêche bien sûr que ces conditions ne soient pas les mêmes sur des parties de frontières différentes. Nous traitons dans ce chapitre le cas habituel où on impose des conditions sur \(\mathit{PRE1}\) et/ou \(\mathit{PRE2}\) , par opposition au chapitre suivant où nous parlerons de relations linéaires entre inconnues.

\({p}_{c}={p}_{\text{gz}}-{p}_{\text{lq}}={p}_{{}_{c}}^{\text{imp}}\)

\({p}_{\text{gz}}={p}_{\text{as}}+{p}_{\text{vp}}={p}_{{}_{\text{gz}}}^{\text{imp}}\)

Les flux sont alors des résultats de calculs par [éq A5.1-3] et [éq A5.1-4]

  • Dirichlet \(\mathit{PRE1}\) , Neuman \(\mathit{PRE2}\)

C’est le cas où on impose une valeur à \(\mathit{PRE1}\) et une valeur au flux associé à \(\mathit{PRE2}\) , en ne disant rien sur \(\mathit{PRE2}\) ou en donnant une valeur à FLUN_HYDR2 de FLUX_THM_REP dans AFFE_CHAR_MECA. Appelons \({M}_{{2}_{\text{ext}}}\) cette quantité imposée, qui vaudra \(0\) si rien n’est dit relativement à \(\mathit{PRE2}\) . Nous noterons \({p}_{1}={p}_{{}_{1}}^{\text{imp}}\) la condition imposée sur \(\mathit{PRE1}\)

Ceci correspond à:

\(\begin{array}{}{p}_{c}={p}_{\text{gz}}-{p}_{\text{lq}}={p}_{{}_{c}}^{\text{imp}}\\ {p}_{{}_{1}}^{\text{imp}}={p}_{{}_{c}}^{\text{imp}}\end{array}\)

Pour faire la démonstration dans le cadre non homogène, il faudrait introduire un relèvement de la condition \({p}_{1}={p}_{{}_{1}}^{\text{imp}}\) (c’est à dire un champ particulier vérifiant cette condition). Cela alourdit les écritures et n’apporte rien, on se place donc dans le cadre homogène \({p}_{{}_{1}}^{\text{imp}}=0\)

Dans [éq A5.1-3] et [éq A5.1-4], on peut donc prendre et \({\pi}_{2}\) quelconque et \({\pi}_{1}\) vérifiant \({\pi}_{1}=0\) sur \(\partial {O}_{d}\) On commence alors à prendre \({\pi}_{1}=0\) et \({\pi}_{2}=0\) sur tout le bord \(\partial \Omega ` et on obtient [:ref:`éq A5.1-1 <éq A5.1-1>\)] et [éq A5.1-2] au sens des distributions. On multiplie alors [éq A5.1-1] par \({\pi}_{1}\) tel que \({\pi}_{1}=0\) sur \(\partial {\Omega}_{d}\) on multiplie [éq A5.1-2] par \({\pi}_{2}\) quelconque, on intègre par partie, on tient compte de [éq A5.1-3] et [éq A5.1-4] et on obtient, en désignant par \(n\) la normale au bord:

\({\int}_{\partial {\Omega}_{d}}{M}_{\text{as}}.n.{\pi}_{2}d\Gamma ={\int}_{\partial {\Omega}_{d}}{M}_{2\text{ext}}{\pi}_{2}d\Gamma ,\phantom{\rule{4em}{0ex}}\forall {\pi}_{2}\)

On en déduit

\(\begin{array}{c}{M}_{\text{as}}.n={M}_{2\text{ext}}\text{sur}\partial {\Omega}_{d}\end{array}\)

  • Dirichlet \(\mathit{PRE}2\) , Neuman \(\mathit{PRE}1\)

C’est le cas où on impose une valeur à \(\mathit{PRE}2\) et une valeur au flux associé à \(\mathit{PRE}1\) , en ne disant rien sur PRE1 ou en donnant une valeur à FLUN_HYDR1 de FLUX_THM_REP dans AFFE_CHAR_MECA. Appelons \({M}_{{1}_{\text{ext}}}\) cette quantité imposée, qui vaudra \(0\) si rien n’est dit relativement à \(\mathit{PRE}2\) . Nous noterons \({p}_{2}={p}_{{}_{2}}^{\text{imp}}\) la condition imposée sur \(\mathit{PRE}2\) .

Ceci correspond à:

\(\begin{array}{c}{p}_{\text{gz}}={p}_{\text{as}}+{p}_{\text{vp}}={p}_{{}_{\text{gz}}}^{\text{imp}}\\ {p}_{2}^{\text{imp}}={p}_{{}_{\text{gz}}}^{\text{imp}}\end{array}\)

La démonstration est la même qu’au paragraphe précédent et aboutit à:

\(\begin{array}{c}({M}_{\text{lq}}+{M}_{\text{vp}}).n={M}_{1\text{ext}}\text{sur}\partial {O}_{d}\end{array}\)

A 2 .3 Cas de conditions aux limites faisant intervenir des relations linéaires entre inconnues principales

Code_Aster permet d’introduire comme conditions aux limites des relations entre degrés de liberté, portés par le même nœud ou par des nœuds différents. Cette possibilité est atteinte via le mot clé LIAISON_DDL de la commande AFFE_CHAR_MECA.

Soit \({p}_{{}_{\text{lq}}}^{\text{imp}}\) la valeur que l’on veut imposer à la pression de liquide sur \(\partial {O}_{d}\) . Compte tenu de [éq A5.1-5], et du choix des inconnues principales pour ce comportement, on écrit:

\({p}_{\text{gz}}-{p}_{c}={p}_{2}-{p}_{1}={p}_{{}_{\text{lq}}}^{\text{imp}}\) éq A5.3-1

Les relations linéaires sont traitées dans Aster par introduction de multiplicateurs de Lagrange. Ceci correspond en l’espèce à la formulation suivante:

Trouver \({p}_{1},{p}_{2,}\mu\) tels que:

\(\begin{array}{}-{\int}_{\Omega}(\dot{{m}_{\text{lq}}}+{\dot{m}}_{\text{vp}}){\pi}_{1}d\Omega +{\int}_{\Omega}({M}_{\text{lq}}+{M}_{\text{vp}}).\nabla {\pi}_{1}d\Omega -{\int}_{\Omega}\dot{{m}_{\text{as}}}{\pi}_{2}d\Omega +{\int}_{\Omega}{M}_{\text{as}}.\nabla {\pi}_{2}d\Omega +\\ +{\int}_{\partial {\Omega}_{d}}{\mu}^{\text{*}}({p}_{2}-{p}_{1}-{p}_{\mathrm{lq}}^{\text{imp}})d\Gamma +{\int}_{\partial {\Omega}_{d}}\mu ({\pi}_{2}-{\pi}_{1}-{p}_{\mathrm{lq}}^{\text{imp}})d\Gamma \forall {\pi}_{1},{\pi}_{2},{\mu}^{\text{*}}\end{array}\) éq A5.3-2

Pour faire la démonstration dans le cadre non homogène, il faudrait introduire un relèvement de la condition \({p}_{2}-{p}_{1}-{p}_{{}_{\text{lq}}}^{\text{imp}}=0\) (C’est à dire des champs particuliers vérifiant cette condition). Cela alourdit les écritures et n’apporte rien, on se place donc dans le cadre homogène \({p}_{{}_{\text{lq}}}^{\text{imp}}=0\) .

On commence alors à prendre \({\pi}_{1}=0\) et \({\pi}_{2}=0\) sur tout le bord \(\partial \Omega\) et on obtient [éq A5.1-1] et [éq A5.1-2] au sens des distributions. On multiplie alors [éq A5.1-1] par \({\pi}_{1}\) quelconque on multiplie [éq A5.1-2] par \({\pi}_{2}\) quelconque, on intègre par partie, on porte les résultats trouvés dans [éq A5.3-2] et on obtient:

\(\begin{array}{}{\int}_{\partial {\Omega}_{d}}({M}_{\text{lq}}+{M}_{\text{vp}}).n.{\pi}_{1}d\Gamma +{\int}_{\partial {\Omega}_{d}}{M}_{\text{as}}.n.{\pi}_{2}d\Gamma +\\ {\int}_{\partial {\Omega}_{d}}{\mu}^{\text{*}}({p}_{2}-{p}_{1})d\Gamma +{\int}_{\partial {\Omega}_{d}}\mu ({\pi}_{2}-{\pi}_{1})d\Gamma =0\forall {\pi}_{1},{\pi}_{2},{\mu}^{\text{*}}\end{array}\) éq A5.3-3

Il est clair que [éq A5.3-3] redonne bien \({p}_{2}-{p}_{1}={p}_{{}_{\text{lq}}}^{\text{imp}}=0\)

En prenant de plus \({\pi}_{2}-{\pi}_{1}=0\) , on trouve:

\({\int}_{\partial {\Omega}_{d}}({M}_{\text{lq}}+{M}_{\text{vp}}+{M}_{\text{as}}).n.{\pi}_{1}d\Gamma =0\forall {\Gamma}_{1}\)

D’où l’on déduit:

\(\begin{array}{c}({M}_{\text{lq}}+{M}_{\text{vp}}+{M}_{\text{as}}).n=0\text{sur}\partial {\Omega}_{d}\end{array}\) éq A5.3-4

A2.4 Cas de conditions aux limites mixtes

Nous abordons dans cette section le cas des conditions aux limites mixtes (ou conditions d’échange) appliquées sous forme d’une condition de Neumann. Ces conditions peuvent porter sur l’une ou l’autre des équations de conservation ou bien sur les deux.

L’écriture des flux se fait alors sous la forme:

\(\begin{array}{c}\left({M}_{\text{lq}}+{M}_{\text{vp}}\right).n=\phantom{\rule{0.5em}{0ex}}{M}_{1\text{ext}}\text{sur}\partial {\Omega}_{d}\end{array}\) avec \({M}_{1\text{ext}}={\lambda}_{11}({p}_{1}-{p}_{{}_{1}}^{\text{ext}})-{\lambda}_{12}({p}_{2}-{p}_{{}_{2}}^{\text{ext}})\)

et/ou

\(\begin{array}{c}{M}_{\text{as}}.n={M}_{2\text{ext}}\text{sur}\partial {\Omega}_{d}\end{array}\) avec \({M}_{2\text{ext}}={\lambda}_{21}({p}_{1}-{p}_{{}_{1}}^{\text{ext}})-{\lambda}_{22}({p}_{2}-{p}_{{}_{2}}^{\text{ext}})\)

Cette condition est implémentée de manière explicite (la valeur de \(\mathit{PRE1}\) et/ou \(\mathit{PRE2}\) étant prise au temps précédent). Ces conditions mixtes sont donc appliquées comme des conditions de Neumann (cf. précédemment).

A 2 . 5 Les cas non linéaires

Nous ne faisons ici qu’aborder des questions plus difficiles consistant à imposer soit la pression de vapeur soit la pression d’air sec. Compte tenu des relations [éq A5.1-7], [éq A5.1-8] et [éq A5.1-9] imposer une valeur sur la pression de vapeur revient à imposer une relation non linéaire sur la pression de liquide; de même qu’imposer une pression d’air sec sur la pression de gaz.

A titre d’exemple, nous abordons le cas d’un pression d’air sec imposée pour un comportement LIQU_VAPE_GAZ, et nous supposons que nous savons écrire la relation non linéaire reliant la pression de vapeur et la pression de liquide.

La relation à imposer est donc:

\({p}_{\text{as}}={p}_{\text{gz}}-{p}_{\text{vp}}={p}_{2}-{p}_{\text{vp}}={p}_{{}_{\text{as}}}^{\text{imp}}\) éq A5.4-1

En différenciant cette relation, on trouvera une condition sur les variations virtuelles de pressions:

\({\text{dp}}_{\text{as}}={\text{dp}}_{\text{gz}}-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{\text{dp}}_{\text{lq}}={\text{dp}}_{\text{gz}}-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}({\text{dp}}_{\text{gz}}-{\text{dp}}_{c})\)

Soit encore

\({\text{dp}}_{\text{as}}={\text{dp}}_{2}-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}({\text{dp}}_{2}-{\text{dp}}_{1})=\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{\text{dp}}_{1}+(1-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}){\text{dp}}_{2}\)

La formulation variationnelle serait alors:

\(\begin{array}{c}-{\int}_{\Omega}\left(\dot{{m}_{\text{lq}}}+{\dot{m}}_{\text{vp}}\right){\pi}_{1}d\Omega +{\int}_{\Omega}\left({M}_{\text{lq}}+{M}_{\text{vp}}\right).\nabla {\pi}_{1}d\Omega -{\int}_{\Omega}\dot{{m}_{\text{as}}}{\pi}_{2}d\Omega +{\int}_{\Omega}{M}_{\text{as}}.\nabla {\pi}_{2}d\Omega +\\ +{\int}_{\partial {\Omega}_{d}}{\mu}^{\text{*}}\left({p}_{2}-{p}_{\text{vp}}-{p}_{\mathit{as}}^{\text{imp}}\right)d\Gamma +{\int}_{\partial {\Omega}_{d}}\mu \left(\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{\pi}_{1}+\left(1-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}\right){\pi}_{2}\right)d\Gamma \phantom{\rule{4em}{0ex}}\forall {\pi}_{1},{\pi}_{2},{\mu}^{\text{*}}\end{array}\) Et on trouverait:

\({\int}_{\partial {\Omega}_{d}}({M}_{\text{lq}}+{M}_{\text{vp}}).n.{\pi}_{1}d\Gamma +{\int}_{\partial {\Omega}_{d}}{M}_{\text{as}}.n.{\pi}_{2}d\Gamma +{\int}_{\partial {\Omega}_{d}}\mu (\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{\pi}_{1}+(1-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}){\pi}_{2})d\Gamma =0\forall {\pi}_{1},{\pi}_{2}\)

En prenant \(\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{\pi}_{1}+(1-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}){\pi}_{2}=0\) on trouverait:

\(\begin{array}{c}\left(1-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}\right)\left({M}_{\text{lq}}+{M}_{\text{vp}}\right).n-\frac{\partial {p}_{\text{vp}}}{\partial {p}_{\text{lq}}}{M}_{\text{as}}.n=0\end{array}\) éq A5.4-2

Approche numérique pour le calcul des paramètres d’entrée de la loi HYDR_TABBAL

Dans cet annexe, on détaille l’approche numérique utilisée pour calculer les paramètres d’entrée, indispensables à l’utilisation de la loi ’HYDR_TABBAL’ (Voir section 2.5). On rappelle que ces paramètres sont: \({S}_{\mathit{BJH}}(\mathit{HR})\) , \({A}_{0}\) , \({\omega}_{\mathit{BJH}}(\mathit{HR})\) , \(t(\mathit{HR})\) et \({\left(\frac{\partial \gamma }{\partial {\epsilon}_{s}}\right)|}_{\mu}\) . Il est conseillé à l’utilisateur d’évaluer ces paramètres selon l’ordre présenté ci-dessous :

  • \({A}_{0}\) est appelée la surface spécifique du matériau par unité de volume (en m²/m³) ouaussi la surface totale des pores par unité de volume.Il existe divers méthodespour estimer la surface spécifique. Parmi lesquelles, les méthodes BET [3]et BJH [4] qui sont des méthodes numériques baséessur des courbesexpérimentales de rétention d’eau. Comme première étape, Il est conseillé à l’utilisateur de calculer la valeur de la surface spécifique \({S}_{\mathit{BET}}\) par la méthode BETselon les équationsA2.1 et A2.2 :

\(\frac{\mathit{HR}}{w\left(\mathit{HR}-1\right)}=\frac{1+\left(C-1\right)\mathit{HR}}{C{w}_{m}}\) éq A2.1

\({S}_{\mathit{BET}}=\frac{{N}_{A}{A}_{m}{V}_{m}}{{V}_{M}}\) éq A2.2

Dans l’équation A2.1, \(w(\mathit{HR})\) est la courbe de rétention d’eau du matériau écrite en terme de teneur en eau. La méthode BET s’applique normalement dans la gamme de \(\mathit{HR}\) situéeentre 5% et 35%. Dans cette gamme, en traçant la courbe \(\frac{\mathit{HR}}{w\left(\mathit{HR}-1\right)}\) fonction de \((\mathit{HR})\) , on obtient une droite. Àpartir de cette dernière, on peut calculer unepente \(\frac{1}{C{w}_{m}}\) et une ordonnéeà l’origine \(\frac{C-1}{C{w}_{m}}\) . Ceci permettrad’évaluer la valeur des paramètres \(C\) et \({w}_{m}\) de l’équation A2.1. Le volume \({V}_{m}\) (en m3/g)de l’équation A2.2, est le volume d’eau correspondant à \({w}_{m}\) . Finalement, \({S}_{\mathit{BET}}\) est calculéeà partir de l’équation A2.2 avec \({N}_{A}\) (=6.022 1023 mol-1)le nombre d’Avogadro, \({A}_{m}\) (=10.6 Å2) la surface occupée par une molécule d’eau et \({V}_{M}\) (=18.10-6 m3/mol) le volume molaire de l’eau.

\({A}_{0}\) est ensuite calculée à partirdu produit de la densité sèche (g/m3) du matériau par la surface spécifique \({S}_{\mathit{BET}}\) (m²/g).

  • \(t(\mathit{HR})\) est l’épaisseur de la couche d’eau adsorbée. Cette épaisseur peut être évaluée par différentesrelations empiriques qu’on peut trouver dans la littérature. Parmi lesquelles, on cite celle de Badmann [2] qui s’écrit comme: \(t\left(\mathit{HR}\right)={K}_{1}+{K}_{2}\bullet \ln(-\ln\left(\mathit{RH}\right))\) avec \({K}_{1}\) et \({K}_{2}\) des coefficients à recalersur une courbe de rétention d’eau du matériau ( \(t(\mathit{HR})={V}_{w}(\mathit{HR})/{S}_{\mathit{BET}}\) et \({V}_{w}(\mathit{HR})\) le volume d’eau calculé à partir de la courbe de rétention). Pour des matériaux cimentaireset une température ambiante, Badmann calcule \({K}_{1}\) = 0,385 nm et \({K}_{2}\) = -0,189 nm.

  • Pour les tableaux \({S}_{\mathit{BJH}}(\mathit{HR})\) et \({\omega}_{\mathit{BJH}}(\mathit{HR})\) , un calcul BJH selon la méthode itérative introduite par Barret,Joyner et Halenda [4] est demandé. Cette méthode est détaillée exhaustivement dans [5].

Comme données d’entrée nécessairespour effectuer le calcul BJH, on cite: la courbe de rétention d’eau du matériau écrite en teneur en eau ou en saturation (fonction de HR), l’épaisseur de la couche adsorbée \(t(\mathit{HR})\) et une valeur de \({A}_{0}\) selon la méthode BET. Il est conseillé avant de procéder par le calcul BJH que la courbe de rétention d’eau soit recalée par un modèle numérique par ex: lemodèle de Van-Genuchten ou autres ou même une simple interpolation linéaire entre les points expérimentaux.

L’utilisateur commence le calcul BJHen estimantle rayon moyendes pores \(\stackrel{´}{{r}_{p,n}}(\mathit{HR})\) qui seront itérativement vidés de leur eau selon le niveau de HR. Ce rayon se calcule selon les équationsA2.3, A2.4 et A2.5avec \({r}_{k}(\mathit{HR})\) calculé selonla loi de Kelvin-Laplace (équation A2.5ou \({\sigma}_{L/G}\) est la tension superficielle de la surface liquide/gas égale à 0,072 N/m).

\(\stackrel{´}{{r}_{p,n}}(\mathit{HR})=\frac{({r}_{p,n}(\mathit{HR})+{r}_{p,n-1}(\mathit{HR}))}{2}\) éq A2.3

\({r}_{p,n}(\mathit{HR})={r}_{k,n}(\mathit{HR})+{t}_{n}(\mathit{HR})\) éq A2.4

\({r}_{k,n}(\mathit{HR})=\frac{-2{\sigma}_{L/G}{M}_{w}^{\mathit{mol}}}{RT{\rho}_{w}\ln(\mathit{HR})}\) éq A2.5

Ensuite, la quantité d’eau désorbée entre deux étapes consécutives n et n-1 est calculé: \(\delta {V}_{w,n}={V}_{w,n}-{V}_{w,n-1}\) avec \({V}_{w,n}\) calculé à partir de la courbe de rétention d’eau pour un niveau de HR.

L’utilisateur passe au calcul du volume \((\delta {V}_{p,n})\) et de la surface \((\delta {a}_{p,n})\) des pores qui seront itérativement vidés avec le niveau de HR. Ceci se fait respectivement par les équations A2.6 et A2.7.

\(\delta {V}_{p,n}={\left(\frac{\stackrel{´}{{r}_{p,n}}}{\stackrel{´}{{r}_{p,n}}-{t}_{n}}\right)}^{2}\bullet [\delta {V}_{w,n}]-\delta {t}_{n}\bullet \sum_{i=1}^{n-1}[\delta {a}_{p,n}\bullet (1-\frac{{t}_{n}}{\stackrel{´}{{r}_{p,n}}})]\) éq A2.6

\(\delta {a}_{p,n}=2\delta {V}_{p,n}/\stackrel{´}{{r}_{p,n}}\) éq A2.7

Finalement, le volume cumulé \(\sum_{i}\delta {V}_{p,i}\) et la surface cumulée \(\sum_{i}\delta {a}_{p,i}\) des pores vidés sont calculés pour chaque niveau de HR. Ces deux grandeurs serviront dans le calcul (équation A2.8 et A2.9) de \({S}_{\mathit{BJH}}(\mathit{HR})\) et \({\omega}_{\mathit{BJH}}(\mathit{HR})\) . Ces deux grandeurs représentent respectivement la fraction volumique des pores remplis d’eau liquide (saturés) et la fraction surfacique des pores non saturés (remplis par du gaz et une couche d’eau adsorbée):

\({S}_{\mathit{BJH}}(\mathit{HR})=1-\sum_{i}\delta {V}_{p,i}/({\varphi}_{0}\bullet {V}_{0})\) éq A2.8

\({\omega}_{\mathit{BJH}}(\mathit{HR})=\sum_{i}\delta {a}_{p,i}/{S}_{\mathit{BET}}\) éq A2.9

Dans l’équation A2.8, \(({\varphi}_{0}\bullet {V}_{0})\) représente le volume total des pores avec \(({\varphi}_{0})\) la porosité du matériau et \(({V}_{0})\) le volume totale de l’échantillon. \({S}_{\mathit{BET}}\) est la surface spécifique du matériau préalablement calculée par BET.

La validité du calcul BJH peut être vérifié, selon [5], en comparant la valeur de \({S}_{\mathit{BET}}\) avec la surface spécifique obtenu par BJH \({S}_{\mathit{tot},\mathit{BJH}}\) (surface cumulée totale des pores vidés). Une différence de moins que 10% est conseilléeentre ces deux grandeurs.

On présente ci-dessous un exemple de script python pour la méthode BJH:

import numpy as np

# ——————————————- #

# Definition des caracteristiques materiaux

# ——————————————- #

poro = 0.082 #Porosite

Vt = (3.14*(9**2)/4)*0.3 #Volume total de l’echantillon en cm3

# ——————————————- #

# Parametres pour la relation deKelvin-Laplace

# ——————————————- #

massemol = 0.018 #Masse molaire de l’eau

R = 8.314 #Constante des gazs parfaits

T = 293.15 #Temperature en K

rho = 1000. #Masse volumique de l’eau

gamma = 235.8 * (1-0.625*(1-T/647.096))*(1-T/647.096)**(1.256) #Tension superficielle en fonction de la temperature

C = (2.* (gamma *0.001) *massemol)*(10**9)/ (R*T*rho)

rhosec = 2.385 #densite seche du materiau

msec = rhosec*Vt #masse seche du materiau en g

Sbet = 27 #surface specifique (en m2/g)prealablement calcule par BET

# ——————————————- #

# Courbe de retention d’eau

# ——————————————- #

path= »/home/…/Inter_lin.txt » #recuperer la courbe experimentale de retention d’eau

HR_meas = tab_col.T[0] #recuperer le tableau des HR testes

HR_meas = HR_meas[::-1]

S_meas = tab_col.T[1] #recuperer le tableau des degrés saturations

S_meas = S_meas[::-1]

w_meas = (S_meas*poro*1)/(rhosec) #calcul de la teneur en eau

n=len(HR_meas)

Facteur_Pc = (rho*R*T/massemol)/(10**6)

Pc = -Facteur_Pc*np.log(HR_meas) #calcul des pressions capillaires selon Kelvin

Pc[0]=0

# ——————————————- #

# Rayon de Kelvin-Laplace

# ——————————————- #

Facteur_rk = (massemol*2*(gamma*0.001))/(rho*R*T)

rk = -Facteur_rk*(10**9)/np.log(HR_meas)

# ——————————————- #

# Epaisseur de la couche adsorbee – Selon Badmann

# ——————————————- #

K1 = 0.385 #en nm

K2 = -0.189 #en nm

t= K1+K2*np.log(-np.log(HR_meas))#epaisseur en nm

# ——————————————- #

# Rayon du pore

# ——————————————- #

r=rk+t #en nm

# ——————————————- #

# Rayon moyen des pores

# ——————————————- #

rm = []

rm.append(r[0])

for i in range(1,n):

rm.append((r[i]+r[i-1])/2)

rm=np.asarray(rm) #en nm

# ——————————————- #

# Contraintes a imposer sur le premier point

# ——————————————- #

#Pres de HR100% plus deau adsorbee/uniquement eau libre

Sl=[] #Fraction de pores satures (remplis d’eau liquide)

Sl.append(S_meas[0])

Sg=[] #Fraction de pores non satures

Sg.append(1-Sl[0])

SumVP=[] #~Volume cumule des pores non satures en cm3

SumVP.append(poro*Vt*Sg[0])

DeltaVP = [] #Volume des pores non saturesen cm3

DeltaVP.append(SumVP[0])

DeltaVL =[] #Volume de liquide desorbe entre deux HR consecutives

nb=DeltaVP[0]/((rm[0]/(rm[0]-t[0]))**2)

DeltaVL.append(nb)

# ——————————————- #

# Volume de liquide desorbe entre deux HR

# ——————————————- #

VL = S_meas * poro *Vt #en cm3

for i in range(1,n):

DeltaVL.append(VL[i-1]-VL[i])

# ——————————————- #

# Premiere etape de desorption

# ——————————————- #

DeltaVP.append(((rm[1]/(rm[1]-t[1]))**2)*DeltaVL[1])#Volumedes pores désaturés ente deux HR en cm3

DeltaAP = [] # surface des pores desaturés ente deux HR en m2

DeltaAP.append((2*DeltaVP[0]*(10**-6))/(rm[0]*(10**-9)))

DeltaAP.append ((2*DeltaVP[1]*(10**-6))/(rm[1]*(10**-9)))

DeltaAP_t=[]

DeltaAP_t.append(0)

DeltaAP_t.append(DeltaAP[0]*(1-t[1]/rm[0]))

# ——————————————- #

# Deuxieme a nieme etape de desorption

# ——————————————- #

for i in range(2,n):

j=1

b=0

while j <= i:

b = b+(DeltaAP[i-j]*(1-t[i]/rm[i-j]))

j=j+1

DeltaAP_t.append(b)

a = ((rm[i]/(rm[i]-t[i]))**2)

deltat=t[i-1]-t[i]

DeltaVP.append(a*(DeltaVL[i]-(deltat*(10**-9)*DeltaAP_t[i]*(10**6))))

DeltaAP.append ((2*DeltaVP[i]*(10**-6))/(rm[i]*(10**-9)))

# ——————————————- #

# Volume cumule des pores non satures en cm3

# ——————————————- #

for i in range(2,n+1):

j=1

c=0

while j<=i:

c = c + DeltaVP[i-j]

j=j+1

SumVP.append(c)

# ——————————————- #

# Surface cumulee des pores non satures en m2

# ——————————————- #

SumAP=[]

SumAP.append(DeltaAP[0])

for i in range(2,n+1):

j=1

c=0

while j<=i:

c = c + DeltaAP[i-j]

j=j+1

SumAP.append(c)

SumAP=np.asarray(SumAP)

# ——————————————- #

# Calcul des fractions SBJH et WBJH

# ——————————————- #

Sl=[] #Fraction de pores satures nommé aussi SBJH

Sl.append(S_meas[0])

Sg=[] #Fraction de pores non satures

Sg.append(1-Sl[0])

for i in range(1,n):

Sg.append(SumVP[i]/(poro*Vt)) #Fraction de pores non satures

Sl.append(1-Sg[i]) #Fraction de pores satures nommée aussi SBJH

co=0

for i in Sl:

if i>0:

co=co+1

SBJHT = SumAP[co-1]/(Vt*(10**(-6))) #Surface specifique totale obtenu par BJH en m2/m3

print « S_BJH = », SBJH, »en m2/m3 »

S = (SBJH * Vt*10**-6)/msec

print « S_BJH = »,S, »en m2/g »

Ecart = (S – Sbet)*100/Sbet

print « Ecart avec BET = « ,Ecart, »% » #Verifier que cest moins que 10%

A0 = (Sbet * msec) /(Vt*10-6)**

print « A0= »,A0, »m2/m3 »

wg=SumAP/(Sbet*msec) #Fraction surfacique des pores non saturés nommée aussi SBJH w_bjh

  • Finalement, \({\left(\frac{\partial \gamma }{\partial {\epsilon}_{s}}\right)|}_{\mu}\) est appelé le terme de Shuttleworth qui est fortement dépendant de l’orientation et de la morphologie de la surface solide. Ce terme est recalé àpartir d’expérience de retrait de séchage. Selon les travaux de G. Tabbal [1],[6], ce terme peut avoir une contribution négligeable au développement des contraintes et donc peut dans certain cas être pris comme nul.