v6.03.160 SSNP160 – Diffusion d’hydrogène dans un acier élastoplastique#
Résumé:
Ce test de mécanique quasi-statique non linéaire simule la diffusion d’hydrogène dans un volume homogène d’acier en fonction des déformations plastiques.
La modélisation A est effectuée avec des éléments plans (modélisation mécanique C_PLAN, modélisation thermique PLAN_DIAG).
La modélisation A est effectuée avec des éléments 3D (modélisation mécanique 3D, modélisation thermique 3D_DIAG).
Solution de référence#
Méthode de calcul utilisée pour la solution de référence#
Dans le document [1], un modèle de diffusion des atomes d’hydrogène dans les aciers est proposé. Il considère deux types de concentration des atomes d’hydrogène :
\({C}_{L}\) est la concentration dans le réseau cristallin (Lattice)
\({C}_{T}\) est la concentration dans les «pièges» ou lacunes (Traps)
Sans reprendre ici toute la démarche suivie par les auteurs de [1], la formulation proposée pour l’équation de diffusion de \({C}_{L}\) est :
\(\frac{{C}_{L}+{C}_{T}(1-{\theta}_{T})}{{C}_{L}}\frac{\partial {C}_{L}}{\partial t}-\nabla \cdot ({D}_{L}\nabla {C}_{L})+\nabla \cdot (\frac{{D}_{L}{C}_{L}{V}_{H}}{RT}\nabla {\sigma}_{H})+{\theta}_{T}\frac{d{N}_{T}}{d{\varepsilon}_{\mathit{eq}}^{p}}\frac{d{\varepsilon}_{\mathit{eq}}^{p}}{\mathit{dt}}=0\) eq1
On constate donc que l’équation de diffusion tient compte du gradient local de la trace des contraintes (contrainte hydrostatique \({\sigma}_{H}=1/3\mathit{tr}(\sigma )\) ) et de la déformation plastique équivalente de Von Mises.
Les relations définissant les différentes quantités sont :
\({\theta}_{L}=\frac{{C}_{L}}{{N}_{L}}\) est le taux d’occupation des sites du réseau cristallin, avec \({N}_{T}\) le nombre de sites par unité de volume. \({N}_{L}\) est un constante estimée à \({N}_{L}=5,1{10}^{29}{m}^{\text{-3}}\) pour le fer \(\alpha\) dans [1].
\({\theta}_{T}=\frac{{C}_{T}}{{N}_{T}}\) est le taux d’occupation des sites des pièges, avec \({N}_{T}\) la densité des pièges, c’est à dire le nombre de sites correspondant à des pièges par unité de volume.
Contrairement à \({N}_{L}\) qui est une constante, \({N}_{T}\) est fonction de la déformation plastique suivant l’expression : \({\log}_{10}({N}_{T})={a}_{1}-{a}_{2}\exp(-{a}_{3}{\varepsilon}_{\mathrm{eq}}^{p})\) , avec : \({a}_{1}=23.26,{a}_{2}=2,33,{a}_{3}=-5,5\) [1].
\({D}_{L}=1,27{10}^{\text{-8}}{m}^{2}/s\) \({V}_{H}=2{10}^{\text{-6}}{m}^{3}\) pour le fer \(\alpha\) à température ambiante, \(R=8,3144J/\mathrm{mol}/K\) est la constante des gaz parfaits, et \(T\) la température en \(°K\) .
Il reste à définir \({C}_{T}\) en fonction de \({C}_{L}\) : d’après [1] à l’équilibre, ce qui est le cas pour la CSC :
\({C}_{T}=\frac{{N}_{T}}{1+\frac{1}{{K}_{T}{\theta}_{L}}}\) , avec \({K}_{T}=\exp(-\Delta {E}_{T}/RT)=4,97{10}^{10}\) à température ambiante, suivant [1], en prenant \(\Delta {E}_{T}=-60\mathrm{KJ}/\mathrm{mol}\) .
De façon similaire à la thermique non linéaire,la formulation variationnelle s’écrit alors:
Soit \(\Omega\) un ouvert de \({R}^{3}\) , de frontière \(\Gamma ={\Gamma}_{1}\cup {\Gamma}_{2}\) .
On doit résoudre l’équation [eq. 1] en \({C}_{L}\) sur \(\Omega \times ]0,t[\) avec les conditions aux limites:
\(\lbrace \begin{array}{ccc}{C}_{L}={C}_{L}^{d}& \text{sur}& {\Gamma}_{1}\\ J.n=\phi & \text{sur}& {\Gamma}_{2}\end{array}\) éq 3
avec \(J=-{D}_{L}\nabla {C}_{L}+\frac{{D}_{L}{V}_{H}}{RT}{C}_{L}\nabla {\sigma}_{H}\)
et avec des conditions initiales \({C}_{L}(t=0)\) (et \({C}_{T}(t=0)\) ).
Soit \(v\) une fonction suffisamment régulière s’annulant sur \({\Gamma}_{1}\) , la formulation variationnelle du problème s’écrit :
\(\begin{array}{c}\underset{\Omega}{\int}{D}^{\text{*}}({C}_{L})\frac{{\mathit{dC}}_{L}}{\text{dt}}vd\Omega +\underset{\Omega}{\int}{D}_{L}\nabla {C}_{L}\cdot \nabla vd\Omega -\underset{\Omega}{\int}\nabla v\cdot (\frac{{D}_{L}{C}_{L}{V}_{H}}{RT}\nabla {\sigma}_{H})d\Omega =\\ \underset{\Omega}{\int}{r}_{\text{vol}}vd\Omega +\underset{{\Gamma}_{2}}{\int}\phi vd{\Gamma}_{2}\end{array}\) éq 4
avec :
\({D}^{\text{*}}({C}_{L})=\frac{{C}_{L}+{C}_{T}(1-{\theta}_{T})}{{C}_{L}}\) \(\phi =({D}_{L}\nabla {C}_{L}-\frac{{D}_{L}{V}_{H}}{RT}{C}_{L}\nabla {\sigma}_{H}).n\) et \({r}_{\mathrm{vol}}=-{\theta}_{T}\frac{d{N}_{T}}{d{\varepsilon}_{\mathrm{eq}}^{p}}\frac{d{\varepsilon}_{\mathrm{eq}}^{p}}{\mathrm{dt}}=0\)
La résolution numérique par éléments finis est donc similaire à celle de la thermique non linéaire, et s’appuie sur une \(\theta\) -méthode, modulo deux particularités :
le terme source \({r}_{\mathrm{vol}}\) est non linéaire, comme dans la modélisation du séchage [R7.01.12], et sera intégré de façon explicite ;
le terme \(\underset{\Omega}{\int}\nabla v\cdot (\frac{{D}_{L}{C}_{L}{V}_{H}}{RT}\nabla {\sigma}_{H})d\Omega\) est non symétrique, et sera lui aussi reporté au second membre, comme dans [1] et intégré de façon explicite.
Cette discrétisation explicite de ces deux termes n’est pas contraignante : en effet la résolution de l’équation [eq. 4] est effectuée à chaque pas de temps, chaînée avec une résolution mécanique.
Par contre le terme \(\underset{\Omega}{\int}\nabla v\cdot (\frac{{D}_{L}{C}_{L}{V}_{H}}{RT}\nabla {\sigma}_{H})d\Omega\) nécessite d’appliquer un gradient de température \(\nabla {T}_{\mathrm{ini}}\) supposé uniforme dans l’élément. Le second membre élémentaire calculé est : \({\int}_{\Omega}\nabla {T}_{\mathrm{ini}}K\nabla vd\Omega\) où \(K\) est le tenseur des conductivités thermiques.
Connaissant les conditions initiales \({C}_{L}^{0}={C}_{L}(t=0)\) et \({C}_{T}^{0}={C}_{T}(t=0)\) , on calcule \({C}_{\mathrm{tot}}^{0}={C}_{L}^{0}+{C}_{T}^{0}\) .
et on adimensionne l’inconnue (concentration); la formulation variationnelle est alors :
\(\begin{array}{c}\underset{\Omega}{\int}{D}^{\text{*}}\frac{{\mathit{dc}}_{L}}{\text{dt}}vd\Omega +\underset{\Omega}{\int}{D}_{L}\nabla {c}_{L}\cdot \nabla vd\Omega =\underset{\Omega}{\int}\nabla v\cdot (\frac{{D}_{L}{V}_{H}}{RT}{c}_{L}\nabla {\sigma}_{H})d\Omega +\underset{\Omega}{\int}{\stackrel{ˉ}{r}}_{\text{vol}}vd\Omega +\underset{{\Gamma}_{2}}{\int}\stackrel{ˉ}{\phi }vd{\Gamma}_{2}\end{array}\)
avec : \({c}_{L}=\frac{{C}_{L}}{{C}_{\mathrm{tot}}^{0}}\) \({\stackrel{ˉ}{r}}_{\mathit{vol}}=\frac{-{\theta}_{T}}{{C}_{\mathit{tot}}^{0}}\frac{d{N}_{T}}{d{\varepsilon}_{\mathit{eq}}^{p}}\frac{d{\varepsilon}_{\mathit{eq}}^{p}}{\mathit{dt}}\) et \(\stackrel{ˉ}{\phi }=(-{D}_{L}\nabla {c}_{L}+\frac{{D}_{L}{V}_{H}}{RT}{c}_{L}\nabla {\sigma}_{H}).n\)
Le coefficient de capacité thermique \({D}^{\text{*}}\) est un champ représenté par une variable de commande (NEUT1) .
Dans un élément de volume isolé (\(\nabla {C}_{L}=0\) sur la frontière ) et chargé de façon à ce que l’état de contraintes et de déformations soit uniforme, la concentration totale \({C}_{\mathrm{tot}}={C}_{l}+{C}_{t}\) est constante.
La solution analytique décrite dans [1] est obtenue directement sous forme de l’équation du second degré : \({C}_{T}^{2}-B{C}_{T}+{N}_{T}{C}_{\mathrm{tot}}=0\) , avec \(B={N}_{L}/{K}_{T}+{C}_{\mathrm{tot}}+{N}_{T}\) : un des deux solutions mènerait à des valeurs négatives de \({C}_{L}\) , la solution est donc . \({C}_{T}=1/2(B-\sqrt{{B}^{2}-4{N}_{T}{C}_{\mathrm{tot}}})\)
La représentation de \({C}_{T}({\varepsilon}_{\mathrm{eq}}^{p})\) et \({C}_{L}({\varepsilon}_{\mathrm{eq}}^{p})\) obtenues analytiquement et numériquement dans [1] est :
Résultats de référence#
Incertitude sur la solution#
Incertitude inférieure à \(\text{1 \%}\) .
Références bibliographiques#
«Hydrogen transport near a blunting crack tip» A.H.M. Krom, R.W.J.Koers, A.Bakker in «Journal of the Mechanics and Physics of Solids» 45 (1999) 971-992
Modélisation A#
Caractéristiques de la modélisation A#
Mécanique: modélisation C_PLAN
Thermique: modélisation PLAN_DIAG
Caractéristiques du maillage#
Un élément QUAD4. Nombre de nœuds : 4
Grandeurs testées et résultats#
La solution analytique étant obtenue en fonction de la déformation plastique, on teste d’abord cette valeur en fonction du temps, puis la concentration en hydrogène en fonction du temps:
Identification |
Instant |
Référence |
Tolérance % |
P1 |
11000000.0 |
0.0198 |
0.10% |
P1 |
11500000.0 |
0.0297 |
0.10% |
P1 |
12000000.0 |
0.0396 |
0.10% |
P1 |
12500000.0 |
0.0495 |
0.10% |
P1 |
13000000.0 |
0.0594 |
0.10% |
P1 |
15100000.0 |
0.100980 |
0.10% |
P1 |
20000000.0 |
0.1980 |
0.10% |
CT2 |
11000000.0 |
0.5058 |
0.10% |
CT2 |
11500000.0 |
0.6522 |
0.3% |
CT2 |
12000000.0 |
0.8239 |
0.4% |
CT2 |
12500000.0 |
0.964 |
1.0% |
CT2 |
13000000.0 |
0.989 |
0.4% |
CT2 |
15100000.0 |
0.99807 |
0.10% |
CT2 |
20000000.0 |
0.999628 |
0.10% |
Modélisation B#
Caractéristiques de la modélisation B#
Mécanique: modélisation 3D
Thermique: modélisation3D_DIAG
Caractéristiques du maillage#
Un élément HEXA8. Nombre de nœuds : 8
Grandeurs testées et résultats#
La solution analytique étant obtenue en fonction de la déformation plastique, on teste d’abord cette valeur en fonction du temps, puis la concentration en hydrogène en fonction du temps:
Identification |
Instant |
Référence |
Tolérance % |
P1 |
11000000.0 |
0.0198 |
0.10% |
P1 |
11500000.0 |
0.0297 |
0.10% |
P1 |
12000000.0 |
0.0396 |
0.10% |
P1 |
12500000.0 |
0.0495 |
0.10% |
P1 |
13000000.0 |
0.0594 |
0.10% |
P1 |
15100000.0 |
0.100980 |
0.10% |
P1 |
20000000.0 |
0.1980 |
0.10% |
CT2 |
11000000.0 |
0.5058 |
0.10% |
CT2 |
11500000.0 |
0.6522 |
0.3% |
CT2 |
12000000.0 |
0.8239 |
0.4% |
CT2 |
12500000.0 |
0.964 |
1.0% |
CT2 |
13000000.0 |
0.989 |
0.4% |
CT2 |
15100000.0 |
0.99807 |
0.10% |
CT2 |
20000000.0 |
0.999628 |
0.10% |
Synthèse des résultats#
La solution calculée est très proche de la solution analytique et valide la méthodologie de simulation de la diffusion d’hydrogène dans un acier.