22 Récipient sphérique sous pression interne
Problème : on étudie le comportement d’un récipient sphérique soumis à une pression interne (un gaz liquéfié) dans l’hypothèse des petites perturbations (HPP). Dans un premier temps, seule la pression interne est appliquée (étude purement élastique).
Données : le récipient est une sphère creuse de rayon intérieur \(R_i = 25\) cm et de rayon extérieur \(R_e = 75\) cm, en matériau élastique linéaire isotrope d’aluminium (\(E = 70\) GPa, \(\nu = 0{,}33\)). Le gaz exerce une pression uniforme \(p_i = 100\) MPa sur la face interne ; la face externe est à la pression atmosphérique de référence \(p_e = 0\).
Paramètres :
| Grandeur | Symbole | Valeur |
|---|---|---|
| Rayon intérieur | \(R_i\) | \(0{,}25\) m |
| Rayon extérieur | \(R_e\) | \(0{,}75\) m |
| Module d’Young | \(E\) | \(70\) GPa |
| Coefficient de Poisson | \(\nu\) | \(0{,}33\) |
| Pression interne | \(p_i\) | \(100\) MPa |
L’objectif est de justifier une modélisation axisymétrique 2D, de mener une étude de convergence du maillage, d’extraire le déplacement maximal, de discuter la forme des tenseurs des contraintes et des déformations aux pôles, et de comparer le déplacement radial le long de l’épaisseur à la solution analytique de Lamé.
22.1 I) Justification du modèle axisymétrique 2D
22.1.1 Géométrie, chargement, matériau
Le récipient est une sphère creuse, la pression interne est uniforme et normale à la surface interne, et le matériau est isotrope. Les trois ingrédients du problème (géométrie, chargement, matériau) sont donc invariants par rotation autour d’un axe quelconque passant par le centre.
22.1.2 Conséquence : champs à symétrie sphérique
Dans un repère sphérique centré \((r, \theta, \varphi)\), toutes les grandeurs sont indépendantes de la latitude \(\theta\) et de la longitude \(\varphi\) : elles ne dépendent que de la coordonnée radiale \(r\). Les directions tangentielles étant équivalentes, on a en tout point
\[\sigma_{\theta\theta} = \sigma_{\varphi\varphi}, \qquad \varepsilon_{\theta\theta} = \varepsilon_{\varphi\varphi},\]
et le seul déplacement non nul est le déplacement radial \(u_r(r)\) : la sphère se dilate uniformément.
22.1.3 Modèle 2D axisymétrique
La symétrie de révolution autour d’un axe (disons l’axe \(z\)) permet de réduire la géométrie à sa section méridienne dans le demi-plan \((r,z)\), \(r \ge 0\) — ici un quart d’anneau compris entre les cercles de rayons \(R_i\) et \(R_e\) (\(r \ge 0\), \(z \ge 0\)). La révolution de cette section engendre la sphère complète. On peut donc répondre oui à la question 1 : le problème se modélise en 2D axisymétrique, pourvu que la formulation intègre correctement la contribution circonférentielle \(\varepsilon_{\theta\theta} = u_r / r\) (effet de l’anneau de révolution).
Remarque — réduction au quart de coupe. La symétrie par rapport au plan équatorial (\(u_z = 0\) sur \(z = 0\)) permet de ne conserver qu’un quart de la section méridienne (\(z \ge 0\)) et de diviser par deux le nombre d’éléments. Le modèle est mathématiquement équivalent (même solution de Lamé). Une première implémentation éléments finis du quart semblait nettement moins précise (convergence \(O(h)\), erreur ~50 fois plus grande), mais l’étude quantitative du notebook recipient_sph_convergence.ipynb a montré qu’il s’agissait d’une condition limite incomplète : la BC \(u_z = 0\) sur la coupe \(z = 0\) n’était imposée qu’aux sommets, laissant libres les milieux d’arête \(\mathbb{P}_2\) de \(u_z\) sur l’équateur artificiel. Avec la BC complète (sommets et milieux), le quart retrouve exactement la convergence et les valeurs du demi-anneau (erreur \(L^2 =
6{,}16\times10^{-4}\) à \(l_c = 0{,}0125\), ordre \(O(h^2)\)). Le notebook numérique utilise donc le quart de section méridienne.
22.2 II) Solution analytique — sphère de Lamé sous pression interne
En symétrie sphérique, l’équilibre local se réduit à
\[\frac{d\sigma_{rr}}{dr} + \frac{2}{r}\left(\sigma_{rr} - \sigma_{\theta\theta}\right) = 0,\]
et les déformations s’écrivent \(\varepsilon_{rr} = du_r/dr\), \(\varepsilon_{\theta\theta} = \varepsilon_{\varphi\varphi} = u_r/r\). La loi de Hooke isotrope donne alors la solution de Lamé pour une sphère creuse sous pression interne \(p_i\) :
\[\boxed{\sigma_{rr}(r) = A - \frac{2B}{r^3}, \qquad \sigma_{\theta\theta}(r) = \sigma_{\varphi\varphi}(r) = A + \frac{B}{r^3}}\]
\[\boxed{u_r(r) = A\,\frac{(1-2\nu)}{E}\,r + B\,\frac{(1+\nu)}{E}\,\frac{1}{r^2}}\]
avec, comme données par l’énoncé :
\[A = p_i\,\frac{R_i^3}{R_e^3 - R_i^3}, \qquad B = \frac{p_i}{2}\,\frac{R_i^3 R_e^3}{R_e^3 - R_i^3}.\]
Les conditions aux limites \(\sigma_{rr}(R_i) = -p_i\) et \(\sigma_{rr}(R_e) = 0\) sont vérifiées exactement : la face interne est en compression radiale (\(-100\) MPa) et la face externe libre.
Déplacement maximal : \(u_r\) est décroissant en \(1/r^2\), le maximum est atteint sur la face interne \(r = R_i\). Le calcul est effectué numériquement ci-dessous.
A = 3.846 MPa
B = 8.113e+05 N.m
sigma_rr(Ri) = -100.00 MPa | sigma_rr(Re) = 0.00 MPa
sigma_tt(Ri) = 55.77 MPa | sigma_tt(Re) = 5.77 MPa
u_r(Ri) = 0.2513 mm | u_r(Re) = 0.0414 mm
sigma_vm max = 155.77 MPa (face interne)
=> Déplacement maximal (analytique) u_max = 0.2513 mm22.3 III) Tenseurs des contraintes et des déformations aux pôles
Aux pôles de la sphère (points de l’axe de révolution, \(r = 0\), \(z = \pm R_e\) ou \(\pm R_i\)), la direction radiale est portée par l’axe \(z\) et les deux directions tangentielles sont équivalentes. En tout point du récipient, et en particulier aux pôles, le tenseur des contraintes est diagonal dans la base sphérique \((e_r, e_\theta, e_\varphi)\) :
\[\boxed{\boldsymbol{\sigma} = \begin{pmatrix} \sigma_{rr} & 0 & 0 \\ 0 & \sigma_{\theta\theta} & 0 \\ 0 & 0 & \sigma_{\varphi\varphi} \end{pmatrix}} \qquad \text{avec} \quad \sigma_{\theta\theta} = \sigma_{\varphi\varphi}.\]
En effet, il n’existe aucune direction tangentielle privilégiée (symétrie sphérique) : aucun cisaillement ne peut apparaître, d’où la forme diagonale, et les deux composantes circonférentielles sont égales. La base sphérique étant locale, la même forme diagonale est valable aux pôles (le tenseur est sphérique, proportionnel en partie à \(\mathbf{I}\)).
De même, le tenseur des déformations est diagonal avec \(\varepsilon_{\theta\theta} = \varepsilon_{\varphi\varphi} = u_r/r\) :
\[\boxed{\boldsymbol{\varepsilon} = \begin{pmatrix} \varepsilon_{rr} & 0 & 0 \\ 0 & \varepsilon_{\theta\theta} & 0 \\ 0 & 0 & \varepsilon_{\theta\theta} \end{pmatrix}}, \qquad \varepsilon_{rr} = \frac{du_r}{dr} = A\,\frac{1-2\nu}{E} - 2B\,\frac{1+\nu}{E}\,\frac{1}{r^3}.\]
22.3.1 Valeurs numériques aux pôles
Aux pôles externes (\(r = R_e\)) : \(\sigma_{rr} = 0\) (paroi libre), \(\sigma_{\theta\theta} = \sigma_{\varphi\varphi} = 5{,}77\) MPa, et le déplacement vaut \(u_r(R_e) \simeq 41{,}4\ \mu\)m. Aux pôles internes (\(r = R_i\)) : \(\sigma_{rr} = -100\) MPa, \(\sigma_{\theta\theta} = \sigma_{\varphi\varphi} = 55{,}77\) MPa.
En notation axisymétrique du plan méridien \((r,z)\), au pôle supérieur la direction radiale physique est \(z\) : on retrouvera donc numériquement \(\sigma_{zz} = \sigma_{rr}\) et \(\sigma_{rr} = \sigma_{\theta\theta}\) (les deux composantes tangentielles dans le plan méridien et hors-plan).
u_max analytique = 0.2513 mm (face interne, r = Ri)
sigma_VM max = 155.77 MPa (face interne)22.4 IV) Résumé des grandeurs clés
| Grandeur | Formule | Valeur |
|---|---|---|
| \(\sigma_{rr}(R_i)\) | \(-p_i\) | \(-100{,}00\) MPa |
| \(\sigma_{rr}(R_e)\) | \(0\) | \(0\) MPa |
| \(\sigma_{\theta\theta}(R_i)\) | \(A + B/R_i^3\) | \(55{,}77\) MPa |
| \(\sigma_{\theta\theta}(R_e)\) | \(A + B/R_e^3\) | \(5{,}77\) MPa |
| \(u_r(R_i) = u_{\max}\) | \(A\frac{1-2\nu}{E}R_i + B\frac{1+\nu}{E R_i^2}\) | \(0{,}2513\) mm |
| \(u_r(R_e)\) | idem en \(R_e\) | \(0{,}0414\) mm |
| \(\sigma_{VM}^{\max}\) | \(\sigma_{\theta\theta}(R_i) - \sigma_{rr}(R_i)\) | \(155{,}77\) MPa |
22.5 Conclusion
La symétrie sphérique du récipient et de son chargement justifie pleinement la modélisation 2D axisymétrique de la section méridienne : tous les champs ne dépendent que de \(r\), les directions tangentielles sont équivalentes et les tenseurs des contraintes et des déformations sont diagonaux en tout point, aux pôles compris. La solution de Lamé fournit le déplacement maximal \(u_{\max} \simeq 0{,}25\) mm sur la face interne et une contrainte équivalente maximale \(155{,}77\) MPa. La partie numérique va vérifier ces valeurs par éléments finis axisymétriques et réaliser l’étude de convergence du maillage demandée.
22.6 II) Application d’une élévation de température uniforme (\(\Delta T = 100\) °C)
À l’état naturel le récipient est en équilibre thermique avec le milieu ambiant, à la température de référence \(T_0 = 25\) °C. On applique maintenant un changement de température constant et uniforme \(\Delta T = T^* - T_0 = 100\) °C. On reprend les questions 1 (justification du modèle axisymétrique), 3 (configuration déformée et déplacement d’une section) et 4 (tenseurs aux pôles), en vérifiant toujours la convergence du maillage.
22.6.1 Justification du modèle axisymétrique (question 1)
L’élévation de température est uniforme et isotrope : la dilatation thermique libre \(\boldsymbol{\varepsilon}^{th} = \alpha\,\Delta T\,\mathbf{I}\) est identique en tout point et n’introduit aucune direction privilégiée. La géométrie sphérique, le chargement thermique et le matériau isotrope restent invariants par rotation autour d’un axe quelconque passant par le centre : le problème reste à symétrie sphérique complète et la modélisation 2D axisymétrique de la section méridienne reste parfaitement justifiée, avec la même formulation (contribution circonférentielle \(\varepsilon_{\theta\theta} = u_r/r\)).
22.6.2 Solution analytique
Dans un corps homogène, isotrope et libre soumis à une élévation uniforme de température, les dilatations thermiques sont compatibles (champ homogène) : il n’apparaît aucune contrainte élastique. La déformation est purement thermique et hydrostatique, le tenseur des contraintes est nul et le déplacement est radial et linéaire en \(r\) (la sphère se dilate uniformément) :
\[\boxed{\boldsymbol{\varepsilon} = \alpha\,\Delta T\,\mathbf{I}}, \qquad \boxed{\boldsymbol{\sigma} = \mathbf{0}}, \qquad \boxed{u_r(r) = \alpha\,\Delta T\, r}\]
22.6.3 Configuration déformée et déplacement d’une section (question 3)
Le déplacement étant proportionnel au rayon, la sphère déformée reste homothétique de la sphère initiale (rapport \(1 + \alpha\,\Delta T\)) : le maillage déformé est la version amplifiée du maillage initial, sans distorsion. Le déplacement d’une section méridienne est maximal sur la face externe :
\[u_r(R_e) = \alpha\,\Delta T\,R_e \simeq 1{,}79\ \text{mm}, \qquad u_r(R_i) = \alpha\,\Delta T\,R_i \simeq 0{,}60\ \text{mm}.\]
22.6.4 Tenseurs aux pôles (question 4)
Partout, et en particulier aux pôles, le tenseur des contraintes est nul et le tenseur des déformations est diagonal et isotrope :
\[\boldsymbol{\sigma} = \mathbf{0}, \qquad \boldsymbol{\varepsilon} = \alpha\,\Delta T \begin{pmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix}.\]
La diagonalité et l’égalité des composantes tangentielles sont triviales ici (aucun cisaillement, aucune direction privilégiée).
22.7 III) Application d’un gradient de température dans l’épaisseur
On impose finalement les températures de paroi à l’équilibre thermique : \(T(R_i) = T_1 = 500\) °C (surface interne au contact du gaz liquéfié) et \(T(R_e) = T_0 = 25\) °C (surface externe au contact de l’air ambiant).
22.7.1 1) Résolution du problème thermique (question 1)
En régime stationnaire, sans source de chaleur, la conduction dans le récipient (conductivité \(k = 237\) W·m⁻¹·°C⁻¹) est régie par l’équation de Laplace. En symétrie sphérique elle se réduit à une équation différentielle du premier ordre dont la solution générale est la forme donnée par l’énoncé :
\[\frac{1}{r^2}\frac{d}{dr}\!\left(r^2\,\frac{dT}{dr}\right) = 0 \quad\Longrightarrow\quad \boxed{T(r) = A + \frac{B}{r}}\]
22.7.2 2) Détermination des constantes A et B (question 2)
Les conditions aux limites \(T(R_i) = T_1\) et \(T(R_e) = T_0\) donnent :
\[B = \frac{T_1 - T_0}{\frac{1}{R_i} - \frac{1}{R_e}}, \qquad A = T_0 - \frac{B}{R_e}.\]
Application numérique (calculée ci-dessous) : \(B \simeq 178{,}13\) °C·m et \(A \simeq -212{,}50\) °C. La température décroît de \(500\) °C sur la face interne à \(25\) °C sur la face externe (gradient radial, indépendant de la conductivité car le champ stationnaire à symétrie sphérique ne dépend que de la géométrie et des températures imposées).
22.7.3 3) Contraintes élastiques induites par le gradient (question 3)
Le champ \(T(r)\) est non uniforme : les dilatations thermiques \(\alpha\,[T(r) - T_0]\) varient avec le rayon et ne sont pas compatibles avec un champ de déplacement. Contrairement au cas de l’élévation uniforme, un corps libre soumis à un gradient de température induit donc bien des contraintes élastiques, même sans action mécanique extérieure. Le problème thermoélastique sphérique (loi de Duhamel–Neumann) se ramène à
\[\frac{d}{dr}\!\left(\frac{1}{r^2}\frac{d}{dr}\!\left(r^2 u_r\right)\right) = \kappa\,T'(r), \qquad \kappa = \frac{1+\nu}{1-\nu}\,\alpha\]
soit, pour \(T(r) = A + B/r\) (intégration directe de \(\int r^2 T\,dr = A r^3/3 + B r^2/2\)) :
\[\boxed{u_r(r) = \kappa\!\left(\frac{A\,r}{3} + \frac{B}{2}\right) + C_1\,r + \frac{C_2}{r^2}}\]
les constantes \(C_1\), \(C_2\) étant fixées par la condition de surface libre \(\sigma_{rr}(R_i) = \sigma_{rr}(R_e) = 0\). Les contraintes radiales et tangentielles s’en déduisent par la loi de Hooke thermoélastique.
22.7.4 Mise en œuvre par éléments finis (question 3)
Ce calcul thermoélastique se réalise en deux étapes découplées dans le notebook numérique FEniCSx :
- une analyse de conduction thermique stationnaire axisymétrique (espace scalaire, températures imposées \(T_1\) et \(T_0\) sur les faces interne et externe) qui fournit le champ \(T(r) = A + B/r\) ;
- une analyse élastique axisymétrique qui réutilise ce champ de température comme terme source thermique dans la loi de Duhamel–Neumann : la contribution \(\gamma\,(T-T_0)\,\mathbf{I}\) est reportée au second membre de la formulation faible.
22.8 IV) Résumé — partie thermique
| Grandeur | Expression | Valeur |
|---|---|---|
| Élévation uniforme | \(\Delta T = T^* - T_0\) | \(100\) °C |
| \(u_r^{\text{uniforme}}(R_i)\) | \(\alpha\,\Delta T\,R_i\) | \(0{,}5975\) mm |
| \(u_r^{\text{uniforme}}(R_e)\) | \(\alpha\,\Delta T\,R_e\) | \(1{,}7925\) mm |
| Tenseurs (uniforme) | \(\boldsymbol{\sigma} = \mathbf{0}\), \(\boldsymbol{\varepsilon} = \alpha\Delta T\,\mathbf{I}\) | — |
| \(A\) (gradient) | \(T_0 - B/R_e\) | \(-212{,}50\) °C |
| \(B\) (gradient) | \((T_1 - T_0)/(1/R_i - 1/R_e)\) | \(178{,}13\) °C·m |
| \(T(R_i)\) / \(T(R_e)\) | \(A + B/r\) | \(500\) / \(25\) °C |
| \(u_r^{\text{grad}}(R_i)\) | formule de Duhamel–Neumann | \(0{,}55\) mm |
| \(u_r^{\text{grad}}(R_e)\) | formule de Duhamel–Neumann | \(1{,}64\) mm |
| \(u_{\max}^{\text{grad}}\) | dans l’épaisseur (vers \(r = 0{,}5\) m) | \(\simeq 1{,}99\) mm |
| \(\sigma_{\theta\theta}^{\text{grad}}(R_i)\) | loi de Hooke thermoélastique | \(\simeq -958\) MPa (compression) |
| \(\sigma_{\theta\theta}^{\text{grad}}(R_e)\) | loi de Hooke thermoélastique | \(\simeq +228\) MPa (traction) |
22.9 Conclusion (partie thermique)
L’élévation uniforme de température ne crée aucune contrainte : la sphère se dilate uniformément (\(u_r \propto r\), déformée homothétique, tenseurs diagonaux isotropes aux pôles). Le gradient radial de température, en revanche, rend les dilatations thermiques incompatibles et induit des contraintes élastiques très élevées (\(\sigma_{\theta\theta}(R_i) \simeq -958\) MPa, soit près d’un gigapascal en compression sur la face interne), bien au-delà de la limite élastique de l’alliage d’aluminium. La partie numérique vérifie ces deux résultats par éléments finis axisymétriques thermomécaniques (FEniCSx) et compare le champ de température à la solution analytique \(T(r) = A + B/r\).
Elevation uniforme de temperature : Delta T = 100 °C (T* = 125 °C)
Deformation thermique : alpha*DeltaT = 2.390e-03
u_r(Ri) = 0.5975 mm | u_r(Re) = 1.7925 mm (max, face externe)
sigma = 0 partout ; eps = alpha*DeltaT*I (diagonal isotrope)
A = -212.50 °C B = 178.12 °C.m
T(Ri) = 500.0 °C T(Re) = 25.0 °C
u_r(Ri) = 0.5458 mm u_r(Re) = 1.6374 mm
sigma_rr : 0 (Ri) / 0 (Re) MPa sigma_tt(Ri) = -958.0 MPa
Contrainte VM maximale = 958 MPa (face interne)