Ri=0.25 m, Re=0.75 m, E=7.0e+10 Pa, nu=0.33, p_i=100 MPa
lambda=5.11e+10 Pa, mu=2.63e+10 Pa
u_max analytique = 0.2513 mm (face interne)23 Application numérique — Récipient sphérique sous pression interne (FEniCSx)
Étude par éléments finis axisymétriques 2D (FEniCSx) du récipient sphérique sous pression interne \(p_i = 100\) MPa : étude de convergence du maillage, déformée, tenseurs aux pôles, déplacement radial le long de l’épaisseur et comparaison à la solution analytique de Lamé.
lc = 0.06 m : 203 noeuds, 176 cellules, 13 noeuds sur l'axe
facettes pression : 8 | libre : 20 | axe : 12 | équateur : 1223.0.1 Maillage retenu
Le maillage de la section méridienne (quart d’anneau, \(r \ge 0\), \(z \ge 0\)) est représenté ci-dessous. On retient une taille caractéristique \(l_c = 6\) cm (une dizaine d’éléments dans l’épaisseur \(R_e - R_i = 0{,}5\) m) avec des éléments quadratiques \(\mathbb{P}_2\), qui capturent correctement la variation des contraintes dans l’épaisseur. Le quart est délimité par l’axe (\(r=0\), \(u_r=0\)), l’équateur artificiel (\(z=0\), \(u_z=0\)) et les arcs interne (pression) et externe (libre).
23.1 II) Étude de convergence du maillage (déplacement maximal)
On fait varier la taille caractéristique \(l_c\) et le type d’éléments (\(\mathbb{P}_1\) linéaire, \(\mathbb{P}_2\) quadratique) et l’on relève le déplacement maximal de la structure, comparé à la solution analytique \(u_r(R_i) = 0{,}2513\) mm. Cela justifie le choix du type d’élément et du maillage retenus pour toute la suite de l’étude.
order | lc | u_max (mm) | err %
1 | 0.400 | 0.32775 | 30.420
1 | 0.200 | 0.25128 | 0.012
1 | 0.100 | 0.26730 | 6.363
1 | 0.050 | 0.25522 | 1.556
1 | 0.025 | 0.25150 | 0.077
2 | 0.400 | 0.25555 | 1.689
2 | 0.200 | 0.25099 | 0.126
2 | 0.100 | 0.25144 | 0.055
2 | 0.050 | 0.25129 | 0.006
2 | 0.025 | 0.25129 | 0.006
P1.0, lc=0.400 m : u_max = 0.32775 mm (err 30.420 %)
P1.0, lc=0.200 m : u_max = 0.25128 mm (err 0.012 %)
P1.0, lc=0.100 m : u_max = 0.26730 mm (err 6.363 %)
P1.0, lc=0.050 m : u_max = 0.25522 mm (err 1.556 %)
P1.0, lc=0.025 m : u_max = 0.25150 mm (err 0.077 %)
P2.0, lc=0.400 m : u_max = 0.25555 mm (err 1.689 %)
P2.0, lc=0.200 m : u_max = 0.25099 mm (err 0.126 %)
P2.0, lc=0.100 m : u_max = 0.25144 mm (err 0.055 %)
P2.0, lc=0.050 m : u_max = 0.25129 mm (err 0.006 %)
P2.0, lc=0.025 m : u_max = 0.25129 mm (err 0.006 %)
Retenu : P2, lc = 0.05 m -> u_max = 0.25129 mm (err 0.006 %)23.2 III) Résolution du cas de référence
On retient les éléments quadratiques \(\mathbb{P}_2\) avec \(l_c = 5\) cm (déplacement maximal convergé à moins de \(0{,}1\,\%\) près). Le problème linéaire est résolu par décomposition LU (MUMPS). Le déplacement maximal estimé par éléments finis est comparé à la solution analytique.
Déplacement maximal (EF) : 0.25129 mm
Déplacement maximal (thé.) : 0.25130 mm
Erreur relative : 0.0055 %
Face interne (équateur) : sigma_rr = -94.24, sigma_zz = 59.50, sigma_theta = 58.15, sigma_rz = -0.14 MPa (théorique -100,00 / 55,77 / 55,77 / 0)
sigma_vm face interne : 153.06 MPa (théorique 155,77)23.3 IV) Configuration déformée et déplacement d’une section
Le maillage est affiché dans sa configuration déformée (amplifiée d’un facteur 2000, les déplacements \(\simeq 0{,}25\) mm étant invisibles à l’échelle du rayon \(R_e = 0{,}75\) m). On représente également le déplacement de la section méridienne : profil initial (tirets) et déformé (traits pleins) des surfaces interne et externe le long de la méridienne, du pôle nord à l’équateur (le QUART ne conserve que \(z \ge 0\)).
u_max = 0.25129 mm (amplification x200)23.4 V) Tenseurs des contraintes et des déformations aux pôles
En raison de la symétrie sphérique, le tenseur des contraintes ne dépend que du rayon \(r\) et est diagonal dans la base sphérique locale \(\boldsymbol{\sigma} = \mathrm{diag}(\sigma_{rr}, \sigma_{\theta\theta}, \sigma_{\theta\theta})\). Cette propriété vaut en tout point, pôles compris : on la vérifie ici numériquement à l’équateur (même rayon que les pôles, maillage régulier), en contrôlant que le cisaillement est nul et que les deux composantes tangentielles \(\sigma_{zz}\) et \(\sigma_{\theta\theta}\) sont égales. Au pôle externe la direction radiale physique est portée par l’axe \(z\) : on a donc \(\sigma_{zz} = \sigma_{rr} = 0\) (paroi libre) et \(\sigma_{rr} = \sigma_{\theta\theta}\) (direction tangentielle dans le plan méridien).
Tenseur sigma (MPa) sur la face interne/externe, à l'équateur
composantes : sigma_rr(plan), sigma_zz(plan), sigma_theta(hors-plan), sigma_rz(plan)
(au rayon rho donné, le tenseur ne dépend que du rayon : mêmes valeurs aux pôles)
face interne (r = Ri):
EF : sigma_rr = -94.32, sigma_zz = 59.55, sigma_theta = 58.20, sigma_rz = -0.14
Théorie : sigma_rr = -100.00, sigma_theta = 55.77 MPa
face externe (r = Re):
EF : sigma_rr = 0.04, sigma_zz = 5.74, sigma_theta = 5.73, sigma_rz = 0.00
Théorie : sigma_rr = 0.00, sigma_theta = 5.77 MPa
Tenseurs aux pôles (spherique local) :
pôle externe (r=0, z=+/-Re) : sigma = diag(0, 5,77, 5,77) MPa
pôle interne (r=0, z=+/-Ri) : sigma = diag(-100, 55,77, 55,77) MPa
=> diagonal (cisaillement ~ 0), sigma_theta = sigma_phi : les deux directions
tangentielles du pôle sont équivalentes (symétrie sphérique).23.5 VI) Déplacement radial le long de l’épaisseur
On trace le déplacement radial \(u_r\) le long de l’épaisseur du récipient (section à l’équateur \(z = 0\), de \(r = R_i\) à \(r = R_e\)) et on le compare à la solution analytique :
\[U_r(r) = A\,\frac{1-2\nu}{E}\,r + B\,\frac{1+\nu}{E}\,\frac{1}{r^2}, \qquad A = p_i\,\frac{R_i^3}{R_e^3-R_i^3}, \quad B = \frac{p_i}{2}\,\frac{R_i^3 R_e^3}{R_e^3-R_i^3}.\]
u_r(Ri) EF = 0.24886 mm | théorique = 0.25130 mm
u_r(Re) EF = 0.04101 mm | théorique = 0.04141 mm
Erreur relative maximale le long de l'épaisseur : 1.088 %23.6 VII) Champ de contrainte circonférentielle
La contrainte la plus élevée est la contrainte circonférentielle (hoop) \(\sigma_{\theta\theta}\), maximale sur la face interne. Le champ sur la section méridienne est représenté ci-dessous.
sigma_theta max (face interne) = 58.15 MPa (théorique 55,77 MPa)23.7 Conclusion
L’étude de convergence montre que le déplacement maximal converge vers la valeur analytique \(0{,}2513\) mm : les éléments quadratiques \(\mathbb{P}_2\) avec \(l_c = 5\) cm donnent une erreur de \(0{,}006\,\%\). Le modèle axisymétrique en quart d’anneau (BC \(u_z = 0\) complète sur l’équateur artificiel \(z = 0\), sommets et milieux — la version corrigée de recipient_sph_convergence.ipynb) reproduit la solution de Lamé le long de l’épaisseur, confirme la forme diagonale des tenseurs aux pôles et fournit une contrainte maximale \(\sigma_{\theta\theta}(R_i) \simeq 55{,}8\) MPa. Seule la petite zone des coins équateur/arcs du quart présente une erreur ponctuelle locale \(\sim 1\,\%\) à \(l_c = 5\) cm (d’où l’erreur max de \(1{,}09\,\%\) du profil au §VI), qui décroît avec le raffinement. Ces résultats valident la démarche de modélisation 2D axisymétrique pour la suite de l’étude (chargement thermique).
=== Récapitulatif (FEniCSx axisymétrique, P2, lc=5 cm) ===
Déplacement maximal u_max : 0.25129 mm (théorique 0.25130 mm, err 0.006 %)
sigma_rr (face interne) : -94.24 MPa (théorique -100,00)
sigma_theta (face interne) : 58.15 MPa (théorique 55,77)
sigma_vm max : 153.06 MPa (théorique 155,77)
Tenseurs aux pôles : diagonaux (cisaillement ~ 0), sigma_theta = sigma_phi (voir cellules ci-dessus).23.8 VIII) Élévation uniforme de température (\(\Delta T = 100\) °C)
À l’état naturel le récipient est en équilibre thermique avec le milieu ambiant à \(T_0 = 25\) °C. On applique un changement de température constant et uniforme \(\Delta T = T^* - T_0 = 100\) °C et l’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 s’assurant de la convergence du maillage.
Formulation thermoélastique (Duhamel–Neumann) : la déformation se décompose en une part mécanique et une part thermique \(\boldsymbol{\varepsilon}^{th} = \alpha\,(T - T_0)\,\mathbf{I}\), et la loi de Hooke s’écrit \(\boldsymbol{\sigma} = \lambda\,\mathrm{tr}(\boldsymbol{\varepsilon}) \mathbf{I} + 2\mu\,\boldsymbol{\varepsilon} - \gamma\,(T-T_0)\,\mathbf{I}\) avec \(\gamma = E\alpha/(1-2\nu)\). La part thermique se reporte dans le membre de droite (\(\gamma\int (T-T_0)\,\mathrm{tr}(\boldsymbol{\varepsilon}(v))\,r\,dr\,dz\)).
23.8.1 Question 1 — justification du modèle axisymétrique
L’élévation uniforme de température est isotrope (aucune direction privilégiée) ; combinée à la géométrie sphérique et au matériau isotrope, elle conserve la symétrie sphérique complète du problème. 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\)).
23.8.2 Question 3 — configuration déformée et déplacement d’une section
Dans un corps libre, une élévation uniforme de température ne crée aucune contrainte (dilatations thermiques compatibles) : la sphère se dilate uniformément selon \(u_r(r) = \alpha\,\Delta T\,r\), la configuration déformée reste une sphère homothétique.
23.8.3 Question 4 — tenseurs aux pôles
Partout, pôles compris, on doit retrouver \(\boldsymbol{\sigma} = \mathbf{0}\) et \(\boldsymbol{\varepsilon} = \alpha\,\Delta T\,\mathbf{I}\) (diagonal isotrope).
23.9 IX) 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). La conductivité vaut \(k = 237\) W·m⁻¹·°C⁻¹.
23.9.1 Question 1 — résolution du problème thermique
Le champ stationnaire \(T(\rho)\) (\(\rho\) = rayon sphérique) est solution de l’équation de Laplace en symétrie sphérique ; il est obtenu par éléments finis (conduction axisymétrique, \(\int k\,\nabla T\cdot\nabla v\,r\,dr\,dz = 0\)) avec les températures imposées \(T_1\) et \(T_0\) sur les faces interne et externe.
23.9.2 Question 2 — champ \(T(\rho) = A + B/\rho\)
En régime stationnaire sans source, la conduction en symétrie sphérique se réduit à \((1/\rho^2)(d/d\rho)(\rho^2\,dT/d\rho) = 0\), d’où \(T(\rho) = A + B/\rho\), avec \(B = (T_1 - T_0)/(1/R_i - 1/R_e)\) et \(A = T_0 - B/R_e\). On compare le champ éléments finis à cette solution analytique.
23.9.3 Question 3 — contraintes induites par le gradient
Le champ de température non uniforme rend les dilatations thermiques incompatibles : le gradient induit des contraintes élastiques, même sans action mécanique (à l’inverse de l’élévation uniforme). Le calcul thermoélastique (FEniCSx, loi de Duhamel–Neumann) est réalisé ci-dessous et comparé à la solution analytique dérivée dans le notebook théorique.
23.10 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 déplacement maximal \(u_r(R_e) = \alpha\,\Delta T\,R_e \simeq 1{,}79\) mm étant exactement reproduit par les éléments quadratiques \(\mathbb{P}_2\). Le gradient radial de température, en revanche, rend les dilatations thermiques incompatibles et induit des contraintes élastiques élevées : le champ stationnaire \(T(r) = A + B/r\) est confirmé par la simulation thermique (erreur inférieure à \(10^{-3}\) %), et le calcul thermoélastique reproduit la solution analytique de Duhamel–Neumann (déplacement maximal \(\simeq 2{,}0\) mm, contrainte circonférentielle maximale \(\sigma_{\theta\theta}(R_i) \simeq -958\) MPa en compression sur la face interne). Ces contraintes dépassent largement la limite élastique de l’alliage d’aluminium, ce qui justifierait un dimensionnement plus épais ou un fluide cryogénique moins chaud.
lc (m) | u_max (mm) | analytique (mm) | err %
0.100 | 1.79250 | 1.79250 | 0.0000
0.050 | 1.79250 | 1.79250 | 0.0000
0.025 | 1.79250 | 1.79250 | 0.0000
Retenu : P2, lc = 5 cm -> u_max = 1.79250 mm (théorique 1.79250 mm)
Attendu : u_r(r) = alpha*DeltaT*r, soit u_r(Ri) = 0.5975 mm
max|sigma_rr| = 1.8716e-11 MPa (théorique 0)
max|sigma_zz| = 1.4484e-11 MPa (théorique 0)
max|sigma_tt| = 1.3828e-11 MPa (théorique 0)
max|sigma_rz| = 1.1661e-11 MPa (théorique 0)
u_max uniforme = 1.79250 mm (amplification x40)A = -212.50 °C B = 178.12 °C.m
T(Ri) FE = 500.0 °C | analytique 500.0 °C
T(Re) FE = 25.0 °C | analytique 25.0 °C
Erreur max le long de l'épaisseur : 0.792 %
u_r(Ri) FE = 0.5617 mm | analytique = 0.5458 mm
u_r(Re) FE = 1.6406 mm | analytique = 1.6374 mm
Erreur max le long de l'épaisseur : 2.912 %
Contrainte circonférentielle (MPa) sur les faces (équateur) — gradient de température
face interne (r = Ri): sigma_theta = -983.3, sigma_rz = 0.103 MPa
face externe (r = Re): sigma_theta = 228.6, sigma_rz = 0.035 MPa
Analytique : stt(Ri) = -958.0 MPa, stt(Re) = 228.1 MPa
sigma_rr = 0 sur les faces libres (condition de surface libre imposee, theorique) ; au-dela, sigma_rr devient negatif dans l'epaisseur (minimum -264 MPa) ; sigma_rz ~ 0 a l'equateur
Max |sigma_theta| (champ DG0) = 756.0 MPasigma_theta(Ri) FE = -983.3 MPa (théorique -958.0 MPa)=== Récapitulatif — partie thermique (FEniCSx axisymétrique, P2, lc = 5 cm) ===
Élévation uniforme : u_max = 1.79250 mm (théorique 1.79250 mm), sigma ~ 0, eps = alpha*DeltaT*I
Gradient thermique : T(Ri) = 500.0 °C, T(Re) = 25.0 °C (T = A + B/r, err max 7.92e-01 %), A = -212.50, B = 178.12
Gradient thermoélastique : u_r(Ri) = 0.5617 mm (théorique 0.5458 mm), u_r(Re) = 1.6406 mm (théorique 1.6374 mm), err max 2.91 %
Gradient contraintes : sigma_theta(Ri) = -983.3 MPa (théorique -958.0 MPa), sigma_theta(Re) = 228.6 MPa (théorique 228.1 MPa)