Modélisation du débit de fuite des enceintes

Couplage des phénomènes THM

Dans le cadre général fortement contraint présenté dans le chapitre précédent, le but de cette thèse est de mettre en place une approche qui privilégie les phénomènes majeurs et qui s’appuie sur une description simplifiée afin de construire un outil d’aide à l’ingénierie. Ce choix conduit à effectuer un certain nombre d’hypothèses pour rester pragmatique. Dans le cadre de l’étude du comportement à long terme des enceintes de confinement, la stratégie de modélisation proposée se concentre sur la phase d’exploitation. L’hypothèse majeure qui en découle est alors l’affaiblissement possible des couplages entre phénomènes pour qu’ils soient traités les uns après les autres (Fig. 2.1). Cette hypothèse a déjà été utilisée dans de nombreux travaux ([49], [93], [58], [33], [85], [53]), et ses principales justifications sont fournies dans cette section. Seul le comportement du béton durci est pris en compte dans ce chapitre (sans phase initiale de comportement au jeune âge). De même, les situations accidentelles avec des variations fortes de pression et température et la présence de vapeur d’eau ne sont pas considérées dans ces travaux.La température est déterminée dans une première étape (section 2.2), indépendamment des autres phénomènes car on néglige les changements d’état ou la circulation forcée de matière dus aux chargements hydriques ou mécaniques (qui restent modérés et dont les vitesses d’évolution sont faibles [33]). Par ailleurs, Granger [49] a démontré que, dans le cadre du séchage des enceintes de confinement, le niveau de précontrainte biaxiale appliqué (12 MPa horizontalement et 8,5 MPa verticalement) et la fissuration de peau du béton due à la dessiccation ont peu d’influence sur l’ensemble du réseau poral de la structure et donc sur ses capacités globales de transfert. C’est alors le champ de saturation qui est calculé dans une deuxième étape (section 2.3) et qui dépend du champ de température précédent. Dans une troisième étape, l’état mécanique de la structure est calculé (section 2.4). Les déformations différées du béton sont divisées en plusieurs composantes, identifiées successivement et, en pratique, considérées indépendantes ([12], [109]). Leurs évolutions sont des fonctions de la température et de la saturation précalculées afin de suivre les constatations expérimentales du Chapitre 1. Enfin, une dernière étape de post-traitement, basée sur les champs de température, de saturation et de déformation désormais disponibles, permet d’obtenir le flux d’écoulement d’air au travers d’un béton à la fois insaturé et fissuré (section 2.5). 

Modèle thermique

 Le transfert de chaleur se fait principalement par conduction thermique dans le volume de béton. Couplée à la loi de Fourier, l’équation utilisée pour calculer la température 𝑇 (en K) est alors la suivante : 𝜌𝑏𝐶𝑝 𝜕𝑇(𝒙,𝑡) 𝜕𝑡 + 𝛁 ∙ [−𝜆𝛁𝑇(𝒙,𝑡)] = 0 (2.1) Où 𝜌𝑏 est la masse volumique du béton (en kg∙m-3 ), 𝐶𝑝 est sa capacité calorifique massique (en J∙kg-1 ∙K -1 ) et 𝜆 est sa conductivité thermique (en W∙m-1 ∙K -1 ). La dépendance de ces paramètres à la température rend l’équation non linéaire. Au bord du domaine considéré, les échanges thermiques avec le milieu environnant se font par convection. −𝜆𝛁𝑇(𝒙,𝑡) = ℎ𝑐 (𝑇(𝒙,𝑡) − 𝑇𝑒𝑥𝑡) (2.2) Où ℎ𝑐 est le coefficient d’échange convectif (en W∙m-2 ∙K -1 ) et 𝑇𝑒𝑥𝑡 la température ambiante du milieu environnant (en K). Le rayonnement peut être modélisé avec une loi similaire à l’équation (2.2) mais n’a pas été introduit dans les applications numériques de ce manuscrit, l’enceinte externe ayant protégé l’enceinte interne au cours de la construction de la maquette. L’équation de diffusion (2.1) est résolue après discrétisation EF grâce au module de thermique non-linéaire disponible dans Code_Aster [37]. 2.3 Modèle hydrique En phase d’exploitation, les conditions d’humidité sont maintenues au-dessus de 40 % dans le volume interne et l’espace entre enceintes d’un BR. Les hypothèses suivantes peuvent alors être considérées pour le béton des EI [36] : – le phénomène de dessiccation est principalement dû aux mouvements d’eau liquide au cœur du béton ; Modèle hydrique 45 Modélisation et prévision du comportement thermo-hydro-mécanique d’une paroi en béton Application au cas des enceintes de confinement des bâtiments réacteurs nucléaires – les cycles d’imbibition/séchage en phase pré-opérationnelle n’affectent que le béton de peau ; – la pression de l’air est négligeable devant la pression liquide. Dans ce cas, le système (1.10) se simplifie en une unique équation de diffusion nonlinéaire (2.3). Le coefficient de diffusion apparent de l’équation dépend de la température précédemment calculée selon une loi de type Arrhenius, comme proposé par Granger [49]. 𝜕𝑆𝑙 (𝒙,𝑡) 𝜕𝑡 + 𝛁 ∙ [ 𝐾𝑖𝑛𝑡 𝑙 𝑘𝑟𝑙(𝑆𝑙 ) 𝜂𝑙𝜙 𝜕𝑃𝑐 (𝑆𝑙 ) 𝜕𝑆𝑙 𝑇(𝒙,𝑡) 𝑇 𝑟𝑒𝑓 𝑒 𝐸𝑎 𝑅 ( 1 𝑇𝑟𝑒𝑓− 1 𝑇 ) 𝛁𝑆𝑙 ] = 0 (2.3) Où 𝑆𝑙 est le degré de saturation en eau (valeurs entre 0 et 1), 𝐾𝑖𝑛𝑡 𝑙 est la perméabilité intrinsèque à l’eau du béton (en m2 ), 𝑘𝑟𝑙 est la perméabilité relative à l’eau (valeurs entre 0 et 1), 𝜂𝑙 est la viscosité dynamique de l’eau liquide (en Pa∙s), 𝜙 est la porosité du béton, 𝑃𝑐 est la pression capillaire (en Pa), 𝑇 𝑟𝑒𝑓 est la température de référence à laquelle toutes les propriétés hydriques sont définies ou mesurées (en K), 𝐸𝑎 est l’énergie d’activation thermique du béton (en J∙mol-1 ), 𝑅 est la constante des gaz parfaits (8,314 J∙mol-1 ∙K -1 ). Le modèle simplifié (2.3) est classiquement utilisé dans la littérature ([49], [117], [58], [33], [36], exprimé aussi en termes de diffusion de l’humidité relative, de la pression capillaire ou de la concentration en eau par exemple) et diffère peu des modèles plus complets tant que les hypothèses initiales sont respectées [36]. Le modèle de rétention de Van Genuchten [116] définit l’évolution de la pression capillaire en fonction de la saturation : 𝑃𝑐 (𝑆𝑙 ) = 𝑃𝑟 (𝑆𝑙 −1⁄𝑚 − 1) 1⁄𝑛 (2.4) Où 𝑃𝑟 (en Pa) et 𝑛 (scalaire supérieur à 1) sont deux paramètres du modèle et 𝑚 = 1 − 1 𝑛 ⁄ . La loi de Mualem [86] décrit ensuite l’évolution de la perméabilité relative à l’eau en fonction de la saturation : 𝑘𝑟𝑙(𝑆𝑙 ) = √𝑆𝑙 [1 − (1 − 𝑆𝑙 1⁄𝑚) 𝑚 ] 2 (2.5) Les conditions aux parois étant plus naturellement exprimées en humidité relative de l’air ambiant, la loi de Kelvin est utilisée pour y relier la pression capillaire : 𝑃𝑐 (ℎ𝑟 , 𝑇) = − 𝜌𝑙𝑅𝑇 𝑀𝑙 ln(ℎ𝑟 ) (2.6) Où 𝜌𝑙 est la masse volumique de l’eau liquide (en kg∙m-3 ), 𝑀𝑙 est la masse molaire de l’eau (en kg∙mol-1 ) et ℎ𝑟 est l’humidité relative (valeurs entre 0 et 1). Le 46 Modélisation du débit de fuite des enceintes Modélisation et prévision du comportement thermo-hydro-mécanique d’une paroi en béton Application au cas des enceintes de confinement des bâtiments réacteurs nucléaires couplage des équations (2.4) et (2.6) permet de retrouver les isothermes de désorption du béton considéré (courbes ℎ𝑟 = 𝑓(𝑆𝑙 ) à température fixée). Grâce à l’analogie entre les équations (2.1) et (2.3), le problème hydrique est également résolu avec le module de thermique non-linéaire disponible dans Code_Aster. 2.4 Modèle mécanique Le modèle proposé considère les principaux phénomènes à l’origine des déformations différées observées dans les structures en béton précontraint. Les câbles de précontrainte sont modélisés par des éléments barre (1D) parfaitement adhérents aux mailles de béton environnantes (3D).

Câbles de précontrainte 

Les pertes de précontrainte dans les câbles sont calculées selon l’ETC-C [1]. Le cadre théorique à disposition dans Code_Aster ([79], [80] et [81]) est utilisé pour modéliser l’évolution de la tension des câbles 𝐹 (en N) en calculant les pertes de précontrainte par recul d’ancrage Δ𝐹𝑎𝑛𝑐, par frottements Δ𝐹𝜇, dues à la relaxation de l’acier Δ𝐹𝑟 et en ajoutant l’effet des déformations instantanée Δ𝐹𝑒𝑙 et différée Δ𝐹𝑑𝑖𝑓 du béton. 𝐹(𝑡, 𝑠) = 𝐹0 − Δ𝐹𝑒𝑙 − Δ𝐹𝑑𝑖𝑓 − Δ𝐹𝜇 − Δ𝐹𝑎𝑛𝑐 − Δ𝐹𝑟 (2.7) Où 𝑡 est le temps (en h), 𝑠 l’abscisse curviligne le long du câble (en m) et 𝐹0 la tension initiale dans le câble (en N). Les pertes instantanées dues à la déformation élastique du béton sont estimées par : Δ𝐹𝑒𝑙(𝑠) = 𝐴𝑠𝐸𝑠Δ𝜎(𝑠) 2𝐸 (2.8) Où 𝐴𝑠 est la section du câble (en m2 ), 𝐸𝑠 le module d’Young de l’acier (en Pa), 𝐸 le module d’Young du béton et Δ𝜎 la contrainte induite dans le béton par la précontrainte (en Pa). Cette dernière est obtenue par une projection des efforts nodaux dans les câbles en efforts internes dans les mailles de béton voisines. Les pertes par frottements sont données par : Δ𝐹𝜇 (𝑠) = 𝐹0[1 − 𝑒 −𝜇(𝛼+𝑘𝑠) ] (2.9) Où 𝜇 est le coefficient de frottement du câble sur le béton (sans unité), 𝛼 la déviation angulaire le long du câble (obtenue, comme l’abscisse curviligne, par une Modèle mécanique 47 Modélisation et prévision du comportement thermo-hydro-mécanique d’une paroi en béton Application au cas des enceintes de confinement des bâtiments réacteurs nucléaires interpolation des coordonnées de la trajectoire du câble) et 𝑘 le coefficient de pertes linéiques (en m-1 ). Le recul d’ancrage Δ modifie le profil de tension 𝐹 ∗ (𝑠) le long du câble dans la zone d’ancrage sur une distance 𝑑 illustrée en Fig. 2.2 et solution du système : { Δ = 1 𝐸𝑠𝐴𝑆 ∫ |𝐹0 − Δ𝐹𝜇 (𝑠) − 𝐹 ∗ (𝑠)|𝑑𝑠 𝑑 0 |(𝐹0 − Δ𝐹𝜇 (𝑠)) 𝐹 ∗ (𝑠)| = |𝐹(𝑑)| 2 (2.10) Fig. 2.2 : Profil de tension avec recul d’ancrage et équilibre pour résoudre le système (2.10) Les pertes dues à la relaxation de l’acier sont calculées par : Δ𝐹𝑟 (𝑡, 𝑠) = 0,8 ∙ 0,66 105 ∙ 𝜌1000𝑒 9,1 𝐹̃(𝑠) 𝐴𝑠𝑓𝑝𝑘 ( 𝑡 1000) 0,75(1− 𝐹̃(𝑠) 𝐴𝑠𝑓𝑝𝑘 ) 𝐹̃(𝑠) (2.11) Où 𝜌1000 est le coefficient de relaxation à 1000 h de l’acier (en %), 𝑓𝑝𝑘 la résistance à rupture de l’acier (en MPa) et 𝐹̃(𝑠) le profil de tension après prise en compte des pertes instantanées, par frottements et par recul d’ancrage. 

Béton 

Le tenseur des déformations totales dans le béton 𝛆(𝒙,𝑡) est décomposé de la manière suivante : 𝛆(𝒙,𝑡) = 𝛆 𝐞𝐥(𝒙,𝑡) + 𝛆 𝐟𝐩(𝒙,𝑡) + 𝛆 𝐟𝐝(𝒙,𝑡) + 𝛆 𝐫𝐝(𝒙,𝑡) + 𝛆 𝐭𝐡(𝒙,𝑡) (2.12) Où 𝛆 𝐞𝐥 est le tenseur des déformations élasto-endommageables, 𝛆 𝐟𝐩 est le tenseur des déformations de fluage propre, 𝛆 𝐟𝐝 le tenseur des déformations de fluage de dessiccation, 𝛆 𝐫𝐝 le tenseur de retrait de dessiccation et 𝛆 𝐭𝐡 le tenseur des dilatations thermiques. Le retrait endogène, ayant lieu au jeune âge, n’est pas pris en compte dans la présente modélisation. 48 Modélisation du débit de fuite des enceintes Modélisation et prévision du comportement thermo-hydro-mécanique d’une paroi en béton Application au cas des enceintes de confinement des bâtiments réacteurs nucléaires La loi de comportement associée à l’équation (2.12) a été développée avec l’outil TFEL/MFront [52], compatible avec Code_Aster comme loi utilisateur externe. Les équations fournies sont exprimées aux points de Gauss de chaque élément fini. 

Dilatation thermique

 La déformation induite par le chargement thermique du béton est considérée isotrope et proportionnelle à la température : 𝛆 𝐭𝐡(𝒙,𝑡) = 𝛼𝑡ℎ(𝑇(𝒙,𝑡) − 𝑇(𝒙, 0))𝐈𝟑 (2.13) Où 𝛼𝑡ℎ est le coefficient de dilatation thermique du béton (en K-1 ) et 𝐈𝟑 est la matrice identité en 3D. 

Retrait de dessiccation De la même manière, la déformation de retrait de dessiccation est isotrope et proportionnelle au degré de saturation du béton [114] : 𝛆 𝐫𝐝(𝒙,𝑡) = 𝜅𝑟𝑑(𝑆𝑙 (𝒙, 0) − 𝑆𝑙 (𝒙,𝑡))𝐈𝟑 (2.14) Où 𝜅𝑟𝑑 est le coefficient de retrait de dessiccation (sans unité).

Fluage propre Comme proposé par Hilaire [53] et validé sur de nombreux résultats expérimentaux ([19], [100], [18] et [53]), le fluage propre est décrit par une chaîne rhéologique de Burger avec patin vieillissant (cf. Fig. 2.3).