v7.31.151 WTNV151 – Prise en compte d’une condition d’échange hydrique en non saturé#

Résumé:

Les tests présentés permettent de vérifier qu’un terme d’échange hydrique au bord est correctement pris en compte dans les modélisations AXIS_*HH2*, D_PLAN_*HH2* et 3D_*HH2* . Sous certaines hypothèses, on réduit le problème à un problème purement hydraulique en pression capillaire. Ce problème s’écrit alors dans le cas linéaire (ECHANGE_THM) comme une équation de la chaleur et peut ainsi être l’analogue d’un problème en thermique. On compare donc la solution du problème hydraulique à la solution du problème thermique, pour qui le terme d’échange est déjà bien pris en compte. Les modélisations d, e, f concernent la validation du terme d’échange non linéaire en humidité relative (ECHANGE_THM_HR).

Solution de référence (modélisations A, B et C)#

Validation du terme d’échange hydrique à partir du terme d’échange thermique#

On pose un certain nombre d’hypothèses afin faire l’analogie avec le cas thermique :

  • la température est constante (pas de calcul thermique);

  • la pression de gaz est constante;

  • La pression de vapeur est négligeables ;

  • On néglige la diffusion de Fick;

  • Le matériau est indéformable (pas de calculs mécanique) ;

  • La gravité est nulle ;

  • L’eau est incompressible \({\rho}_{w}=\mathit{cste}\) .

On rappelle que la loi de Darcy (voir r7.01.11) pour la phase gazeuse s’écrit:

(4939)#\[{M}_{\mathit{gz}}=-{\rho}_{\mathit{gz}}{\lambda}_{\mathit{gz}}^{H}\nabla {p}_{\mathit{gz}}\]

avec:

  • \({M}_{\mathit{gz}}\) le flux degaz;

  • \({\rho}_{\mathit{gz}}\) la masse volumique du gaz;

  • \({\lambda}_{\mathit{gz}}^{H}\) la conductivité hydraulique du gaz;

  • \({p}_{\mathit{gz}}\) la pression de gaz.

Or:

(4940)#\[{p}_{\mathit{gz}}=\mathit{cste}\Rightarrow \nabla {p}_{\mathit{gz}}=0\]

Donc:

(4941)#\[{M}_{\mathit{gz}}={M}_{\mathit{as}}+{M}_{\mathit{vp}}=0\]

Comme on néglige la diffusion de Fick, alors:

(4942)#\[{M}_{\mathit{ad}}={M}_{\mathit{vp}}=0\]

Par (5004) et (5003) :

(4943)#\[{M}_{\mathit{as}}=0\]

Ainsi on peut simplifier le problème qui se résume à l’équation de conservation de la masse d’eau liquide:

(4944)#\[\dot{{m}_{w}}+\text{Div}({M}_{w})=0\]

Or \({m}_{w}={\rho}_{w}\varphi S(1+\mathit{Tr}(\epsilon ))-{\rho}_{w}^{0}{\varphi}^{0}{S}^{0}\) , avec:

  • \({\rho}_{w}\) la masse volumique de l’eau, \({\rho}_{w}^{0}\) la masse volumique à l’état initial ;

  • \(\varphi =\frac{{V}_{\mathit{vide}}}{V}\) la porosité du matériau, \({V}_{\mathit{vide}}\) correspondant aux pores non remplis, \(V\) à la porosité du matériau à l’état initial ;

  • \(S=\frac{{V}_{\mathit{lq}}}{{V}_{\mathit{vide}}}\) la saturation en eau, \({S}^{0}\) la saturation initiale;

  • \(\epsilon\) la déformation du matériau;

On considère que le matériau est indéformable: \(\epsilon =0\) et \(\varphi =\mathit{cste}={\varphi}^{0}\) et que l’eau est incompressible, \({\rho}_{w}={\rho}_{w}^{0}\) . Alors :

(4945)#\[{m}_{w}(t)={\rho}_{w}^{0}\varphi (S(t)-{S}^{0})\]

Ce qui implique:

(4946)#\[\dot{{m}_{w}}(t)={\rho}_{w}\varphi {\partial}_{t}S({p}_{c}(t))={\rho}_{w}\varphi {S}^{'}({p}_{c})\times \dot{{p}_{c}}(t)\]

Comme on néglige les effets de la gravité, la loi de Darcy pour la phase liquide s’écrit:

(4947)#\[{M}_{\mathit{lq}}={M}_{w}=-{\rho}_{w}{\lambda}_{\mathit{lq}}^{H}(S)\nabla {p}_{w}\]

Or \({p}_{c}={p}_{\mathit{gz}}-{p}_{\mathit{lq}}={p}_{\mathit{atm}}-{p}_{w}\) . Alors \(\nabla {p}_{c}=-\nabla {p}_{w}\) et donc \({M}_{w}={\rho}_{w}{\lambda}_{w}^{H}(S)\nabla {p}_{c}\) .

En remplaçant \({M}_{w}\) et \(\dot{{m}_{w}}\) par leurs dernières expressions dans () et en supposant \({\lambda}_{w}^{H}(S)=\frac{K\times {K}_{\mathit{rl}}(S)}{\mu}\) constante:

(4948)#\[{\rho}_{w}\varphi {S}^{'}({p}_{c})\times \dot{{p}_{c}}(t)+{\lambda}_{w}^{H}(S)\text{Div}({\rho}_{w}.\nabla {p}_{c})=0\]

L’eau est supposée incompressible donc:

(4949)#\[{\rho}_{w}\varphi {S}^{'}({p}_{c})\times \dot{{p}_{c}}(t)+{\rho}_{w}{\lambda}_{w}^{H}(S)\text{Div}(\nabla {p}_{c})=0\]

\({\rho}_{w}\varphi {S}^{'}({p}_{c})\times \dot{{p}_{c}}(t)+{\rho}_{w}{\lambda}_{w}^{H}(S)\text{Div}(\text{}\nabla {p}_{c})=0\)

Avec l’opérateur laplacien:

(4950)#\[{\rho}_{w}\varphi {S}^{'}({p}_{c})\times \dot{{p}_{c}}(t)+{\rho}_{w}{\lambda}_{w}^{H}(S)\Delta {p}_{c}=0\]

Finalement:

(4951)#\[{\partial}_{t}{p}_{c}+\frac{{\lambda}_{w}^{H}(S)}{\varphi {S}^{'}({p}_{c})}\Delta {p}_{c}=0\]

Pour fermer le système, on ajoute:

  • lacondition initiale : \({p}_{c}(x,0):=10000\mathit{Pa}\) ;

  • lacondition mixte sur [BC] (ou [BCFG] selon la géométrie) : \({\lambda}^{H}\nabla ({p}_{c}.n)={h}^{H}({p}_{\mathit{ext}}-{p}_{c})\) .

Cette formulation s’apparente à une équation de la chaleur sans terme source, de la forme:

(4952)#\[\rho {C}_{p}{\partial}_{t}T+\text{Div}(q)=0\]

En prenant la loi de Fourier qui relie le flux de chaleur \(q\) au gradient de température, à savoir:

(4953)#\[q=-\lambda \nabla T\]

En supposant que \(\lambda\) est une constante, caractéristique du matériau:

(4954)#\[{\partial}_{t}T-\frac{\lambda}{\rho {C}_{p}}\Delta T=0\]

Où \(\rho\) est la masse volumique du béton, \({C}_{p}\) la chaleur spécifique du matériau et \(\lambda\) la conductivité thermique.

Le problème thermique est ainsi déjà traité dans THER_NON_LINE. On peut alors comparer les résultats du problème thermique avec ceux du problème THM, en prenant des conditions de fermeture de même type que dans le cas THM i.e.en posant :

\(\begin{array}{c}T(x,0)=\mathit{cste}\\ \lambda (\text{}\nabla T.n)={h}_{T}({T}_{\mathit{ext}}-T)\mathit{au}\mathit{bord}\end{array}\)

En s’assurant de plus que:

  • \(T(x,0)=\mathit{cste}:={p}_{c}(x,0)=10000\) ;

  • \(-\rho {C}_{p}≔\varphi {S}^{'}({p}_{c})\)

  • \(\lambda :={\lambda}^{H}\)

  • \({h}_{T}:={h}^{H}\)

  • \({T}_{\mathit{ext}}:={p}_{\mathit{ext}}\)

Par la suite, on s’appuie alors sur les modélisations thermiques AXIS, PLAN_DIAG, et 3D que l’on compare aux modélisations AXIS_THH2MS, D_PLAN_THH2MS et 3D_THH2MS, respectivement.

Modélisation A#

Caractéristiques de la modélisation A#

Les résultats présentés sont issus de la modélisation en 3D_THH2MS. Ils sont comparés à la solution obtenue en thermique linéaire.

Pour compléter, les calculs sont également réalisés en 3D_HH2MS, en 3D_THH2S et en en 3D_HH2S.

Les déplacements et températures étant bloqués, les résultats doivent être les mêmes dans tous les cas. Le but est ici de vérifier la bonne prise en compte du terme d’échange.

Le maillage est composé de 640 éléments HEXA20 .

Valeurs testées et résultats#

On présente des profils de pression capillaire et de température pour trois instants et on vérifie bien que les résultats sont identiques.

../../../../_images/100000000000034A00000253D77848FCCB23B4C0.png

Fig. 768 Comparaison THM/thermique en 3D avec un terme d’échange hydrique#

On effectue des tests sur trois valeurs:

Points \((x,y)\)

Temps ( \(s\) )

PRE1

\((0,2;0)\)

60 s

13744

1800 s

31758

3600 s

41105

Pour des raisons de gain de temps les modélisations autres que 3D_THH2MS ne seront testées que pour 60 s.

Modélisation B#

Caractéristiques de la modélisation B#

Il s’agit du même calcul que précédemment mais avec une modélisation D_PLAN. Le maillage est composé de 50 éléments QUAD8 et 102 éléments SEG3.

Valeurs testées et résultats#

On présente des profils de pression capillaire et de températurepour trois instants et on vérifie bien que les résultats sont identiques.

../../../../_images/100000000000034A000002538FC7A4480A080936.png

Fig. 769 Comparaison THM/thermique en déformations planesavec un terme d’échange hydrique#

Les résultats en * HH2* ont pour référence les résultats du calcul thermique. On vérifie les valeurs du bord droit du barreau sur lequel la condition d’échange est appliquée.

On effectue des tests sur trois valeurs:

Points \((x,y)\)

Temps (\(s\) )

PRE1

\((0,2;0)\)

60 s

13455

1800 s

31919

3600 s

41087

Pour des raisons de gain de temps les modélisations autres que D_PLAN_THH2MS ne seront testées que pour 60 s.

Modélisation C#

Caractéristiques de la modélisation C#

Il s’agit du même calcul que précédemment mais avec une modélisation axisymétrique.

Valeurs testées et résultats#

On présente des profils de pression capillaire et de températurepour 3 instants et on vérifie bien que les résultats sont identiques.

../../../../_images/100000000000034A000002538823F805066C2234.png

Fig. 770 Comparaison THM/thermique en axis-symétrie avec un terme d’échange hydrique#

Les résultats en HH2 ont pour référence les résultats du calcul thermique. On vérifie les valeurs du bord droit du barreau sur lequel la condition d’échange est appliquée.

On effectue des tests sur trois valeurs:

Points \((x,y)\)

Temps (\(s\) )

PRE1

\((0,2;0)\)

60 s

13888

1800 s

32831

3600 s

42829

Pour des raisons de gain de temps les modélisations autres que AXIS_THH2MS ne seront testées que pour \(60s\) .

Modélisation D#

Caractéristiques de la modélisation D#

Les résultats présentés sont issus de la modélisation en 3D_THH2MS. Il s’agit de la même modélisation que la modélisation A mais avec un terme d’échange non linéaire en HR et non en pression.

Pour compléter, les calculs sont également réalisés en 3D_HH2MS, en 3D_THH2S et en en 3D_HH2S.

Les déplacements et températures étant bloqués, les résultats doivent être les mêmes dans tous les cas. Le but est ici de vérifier la bonne prise en compte de ce terme d’échange.

Le maillage est composé de 640 éléments HEXA20 .

Valeurs testées et résultats#

On vérifie que le flux d’eau (“FH11X”+”FH12X” récupéré par SIEF_NOEU ) obtenu sur le bord d’imposition correspond bien à son expression analytique, à savoir \(q_{w}^{ext}{.}n= \alpha.\rho_{vp}.\left( exp\left[ \frac{Pc.M_{H2O}}{\rho_{l}.R.T}\right] -{HR}^{ext} \right)\). Une tolérance de 10% est permise du fait de la taille du maillage et de l’erreur générée par l’interpolation au noeud.

Points \((x,y)\)

Temps ( \(s\) )

“FH11X”+”FH12X”

Valeur analytique

\((0,2;0)\)

1800

7.3754E-09

8.09e-09

Pour des raisons de gain de temps les modélisations autres que 3D_THH2MS ne seront testées que pour 60 s. On vérifie que les pressions à 60 s sont toutes les mêmes quelles que soient les modélisations.

Points \((x,y)\)

Temps ( \(s\) )

PRE1

\((0,2;0)\)

1800

10017

Modélisation E#

Caractéristiques de la modélisation E#

Il s’agit de la même modélisation que la modélisation D mais en D_PLAN. Le maillage est composé de 50 éléments QUAD8 et 102 éléments SEG3. Sont testées successivement les modélisations D_PLAN_THH2MS, D_PLAN_HH2MS, D_PLAN_THH2S et D_PLAN_HH2S.

Valeurs testées et résultats#

Comme précédemment, on vérifie que le flux d’eau (“FH11X”+”FH12X” récupéré par SIEF_NOEU ) obtenu sur le bord d’imposition correspond bien à son expression analytique, à savoir \(q_{w}^{ext}{.}n= \alpha.\rho_{vp}.\left( exp\left[ \frac{Pc.M_{H2O}}{\rho_{l}.R.T}\right] -{HR}^{ext} \right)\). Une tolérance de 10% est permise du fait de la taille du maillage et de l’erreur générée par l’interpolation au noeud.

Points \((x,y)\)

Temps ( \(s\) )

“FH11X”+”FH12X”

Valeur analytique

\((0,2;0)\)

1800

7.3754E-09

8.09e-09

Pour des raisons de gain de temps les modélisations autres que D_PLAN_THH2MS ne seront testées que pour 60 s. On vérifie que les pressions à 60 s sont toutes les mêmes quelles que soient les modélisations.

Points \((x,y)\)

Temps ( \(s\) )

PRE1

\((0,2;0)\)

1800

10019

Modélisation F#

Caractéristiques de la modélisation F#

Il s’agit de la même modélisation et du même maillage que la modélisation E mais en AXI. Sont testées successivement les modélisations AXI_THH2MS, AXI_HH2MS, AXI_THH2S et AXI_HH2S.

Valeurs testées et résultats#

Comme précédemment, on vérifie que le flux d’eau (“FH11X”+”FH12X” récupéré par SIEF_NOEU ) obtenu sur le bord d’imposition correspond bien à son expression analytique, à savoir \(q_{w}^{ext}{.}n= \alpha.\rho_{vp}.\left( exp\left[ \frac{Pc.M_{H2O}}{\rho_{l}.R.T}\right] -{HR}^{ext} \right)\). Une tolérance de 10% est permise du fait de la taille du maillage et de l’erreur générée par l’interpolation au noeud.

Points \((x,y)\)

Temps ( \(s\) )

“FH11X”+”FH12X”

Valeur analytique

\((0,2;0)\)

1800

7.3754E-09

8.09e-09

Pour des raisons de gain de temps les modélisations autres que AXI_THH2MS ne seront testées que pour 60 s. On vérifie que les pressions à 60 s sont toutes les mêmes quelles que soient les modélisations.

Points \((x,y)\)

Temps ( \(s\) )

PRE1

\((0,2;0)\)

1800

10019

Synthèse des résultats#

Dans le cas d’un échange en pression, les solutionsissues des modélisations AXIS_*HH2* , D_PLAN_*HH2* et 3D_*HH2* sont semblables aux solutionsde référence issues des modélisations AXIS, PLAN_DIAG, et 3D, respectivement. Le terme d’échange en hydraulique est donc correctement représenté. Dans le cas non linéaire d’un échange en humidité relative (ou densité), on vérifie que le flux d’eau retrouvé correspond bien au flux analytique.