12  Cuve sphérique sous pression et gradient thermique

Réservoir sphérique pour le stockage de gaz liquéfié — analyse thermoélastique.

On étudie l’état de contrainte et de déplacement dans un réservoir sphérique de rayon intérieur \(R_i = 1\) m et d’épaisseur \(h = 10\) cm, soumis à une pression interne \(p_i = 350\) bar et à un gradient de température \(\Delta T = 100\) °C dans l’épaisseur. La symétrie sphérique du problème permet une modélisation 1D (coordonnée radiale \(r\)).

12.1 Géométrie et symétries

Le réservoir est une sphère creuse de rayon intérieur \(R_i\) et extérieur \(R_e = R_i + h\). Le problème présente une symétrie sphérique complète : - toutes les grandeurs ne dépendent que de \(r\) (coordonnée radiale) ; - les champs de contrainte et déformation sont isotropes dans les directions tangentielles : \(\sigma_{\theta\theta} = \sigma_{\phi\phi}\), \(\varepsilon_{\theta\theta} = \varepsilon_{\phi\phi}\) ; - la modélisation axisymétrique 2D (méridienne) est suffisante en éléments finis ; - on utilise la symétrie de révolution autour de l’axe \(z\).

12.1.1 Paramètres

Grandeur	Symbole	Valeur
Rayon intérieur	\(R_i\)	\(1.0\) m
Épaisseur	\(h\)	\(0.1\) m
Pression interne	\(p_i\)	\(350\) bar (\(35\) MPa)
Module d’Young	\(E\)	\(200\) GPa
Coefficient de Poisson	\(\nu\)	\(0.3\)
Dilatation thermique	\(\alpha\)	\(1.2\times10^{-5}\) °C\(^{-1}\)
Conductivité thermique	\(k\)	\(50.2\) W/m\(\cdot\)K
Température interne	\(T_i\)	\(-50\) °C
Température externe	\(T_e\)	\(50\) °C
Contrainte admissible	\(\sigma_{\max}\)	\(200\) MPa

12.2 Schéma de la structure

La figure ci-dessous représente la section méridienne de la cuve sphérique avec les conditions aux limites et les chargements.

import matplotlib.pyplot as plt
import numpy as np

def plot_cuve(Ri=1.0, h=0.1):
    Re = Ri + h
    bg = (40/255, 43/255, 48/255, 1.0)
    fig, ax = plt.subplots(figsize=(7, 7))
    fig.patch.set_facecolor(bg); ax.set_facecolor(bg)

    # Demi-cercle intérieur et extérieur
    theta = np.linspace(0, np.pi, 200)
    ax.plot(Ri*np.cos(theta), Ri*np.sin(theta), 'grey', linewidth=2)
    ax.plot(Re*np.cos(theta), Re*np.sin(theta), 'grey', linewidth=2)

    # Hachurage de l'épaisseur
    for t in np.linspace(0.1, 2.9, 12):
        r = Ri + h/3
        ax.plot([Ri*np.cos(t), Re*np.cos(t)], [Ri*np.sin(t), Re*np.sin(t)],
                color='lightgrey', linewidth=0.5, alpha=0.4)

    # Axes
    ax.arrow(0, 0, 1.4*Re, 0, head_width=0.03, head_length=0.04, color='grey')
    ax.arrow(0, 0, 0, 1.4*Re, head_width=0.03, head_length=0.04, color='grey')
    ax.text(1.45*Re, -0.05, "$z$", color='grey', fontsize=14)
    ax.text(-0.08, 1.45*Re, "$r$", color='grey', fontsize=14)

    # Cotes
    ax.annotate("", xy=(0, -Ri), xytext=(0, -Re),
                arrowprops=dict(arrowstyle="<->", color='grey', lw=1.5))
    ax.text(-0.15, -(Ri+Re)/2, "$h$", color='grey', fontsize=13, rotation=90)
    ax.annotate("", xy=(0, 0), xytext=(0, Ri),
                arrowprops=dict(arrowstyle="<->", color='dodgerblue', lw=1.5))
    ax.text(-0.15, Ri/2, "$R_i$", color='dodgerblue', fontsize=13, rotation=90)

    # Pression interne
    for angle in [0.2, 0.5, 0.8, 1.1, 1.4, 1.7, 2.0, 2.3, 2.6]:
        x = (Ri - 0.08)*np.cos(angle)
        y = (Ri - 0.08)*np.sin(angle)
        ax.arrow(x, y, 0.06*np.cos(angle), 0.06*np.sin(angle),
                 head_width=0.03, head_length=0.03, color='dodgerblue')
    ax.text(-0.05, -1.2, r"$p_i$", color='dodgerblue', fontsize=14)

    # Conditions aux limites
    ax.plot([-Re, -Re+0.03], [0, 0.03], 'lightgreen', linewidth=2)
    ax.plot([Re-0.03, Re], [0, 0], 'lightgreen', linewidth=2)
    ax.text(-1.2*Re, 0.05, r"$u_r = 0$", color='lightgreen', fontsize=11)

    ax.set_xlim(-1.6*Re, 1.6*Re)
    ax.set_ylim(-1.6*Re, 1.6*Re)
    ax.set_aspect("equal"); ax.axis("off")
    plt.tight_layout(); plt.show()

plot_cuve()

12.3 Champ de température

L’équation de la chaleur en régime stationnaire dans un milieu isotrope à conductivité constante se réduit au Laplacien nul. En coordonnées sphériques avec symétrie radiale :

\[\nabla^2 T = \frac{1}{r^2}\frac{d}{dr}\left(r^2\frac{dT}{dr}\right) = 0\]

Solution :

\[\boxed{T(r) = A + \frac{B}{r}}\]

Avec les conditions \(T(R_i) = T_i\) et \(T(R_e) = T_e\) :

\[\boxed{T(r) = T_i + (T_e - T_i)\,\frac{1/R_i - 1/r}{1/R_i - 1/R_e}}\]

Pour une paroi mince (\(h \ll R_i\)), le champ est quasi-linéaire dans l’épaisseur.

12.4 Thermoélasticité linéaire

La loi de comportement de Duhamel-Neumann couple les champs mécanique et thermique :

\[\boldsymbol{\sigma} = \mathbb{C} : (\boldsymbol{\varepsilon} - \alpha\,\Delta T\,\mathbf{I})\]

En symétrie sphérique, les déformations s’écrivent :

\[\varepsilon_{rr} = \frac{du_r}{dr}, \qquad \varepsilon_{\theta\theta} = \varepsilon_{\phi\phi} = \frac{u_r}{r}\]

L’équation d’équilibre en l’absence de forces volumiques :

\[\frac{d\sigma_{rr}}{dr} + \frac{2}{r}(\sigma_{rr} - \sigma_{\theta\theta}) = 0\]

En exprimant \(\sigma_{rr}\) et \(\sigma_{\theta\theta}\) en fonction de \(u_r(r)\) via la loi de Hooke thermoélastique, on obtient une équation différentielle pour \(u_r(r)\) dont la solution générale est :

\[u_r(r) = C_1 r + \frac{C_2}{r^2} + \frac{\alpha}{r^2}\int_{R_i}^{r} T(\xi)\,\xi^2\,d\xi\]

12.4.1 Cas (a) — Élévation uniforme de température

Si \(T(r) = T_{\text{amb}}\) constante (dilatation libre sans gradient) : - \(\boldsymbol{\sigma} = \mathbf{0}\) (aucune contrainte) ; - \(u_r(r) = \alpha\,T_{\text{amb}}\,r\) (dilatation libre).

12.4.2 Cas (b) — Gradient thermique \(T(r)\)

Le champ de température non-uniforme induit des contraintes thermoélastiques. Les solutions analytiques complètes pour \(\sigma_{rr}(r)\), \(\sigma_{\theta\theta}(r)\), \(u_r(r)\) sont obtenues par intégration de l’équation différentielle ci-dessus.

12.5 Pression interne seule (sans thermique)

Sous pression interne \(p_i\) seule, la solution de Lamé pour une sphère creuse s’écrit :

\[\boxed{\sigma_{rr}(r) = \frac{p_i R_i^3}{R_e^3 - R_i^3}\left(1 - \frac{R_e^3}{r^3}\right)}\]

\[\boxed{\sigma_{\theta\theta}(r) = \frac{p_i R_i^3}{R_e^3 - R_i^3}\left(1 + \frac{R_e^3}{2r^3}\right)}\]

\[\boxed{u_r(r) = \frac{p_i R_i^3}{R_e^3 - R_i^3}\,\frac{r}{E}\left[(1-2\nu) + (1+\nu)\frac{R_e^3}{2r^3}\right]}\]

12.5.1 Contrainte de Von Mises

Pour un état de contrainte sphérique (\(\sigma_{\phi\phi} = \sigma_{\theta\theta}\)) :

\[\boxed{\sigma_{VM} = |\sigma_{\theta\theta} - \sigma_{rr}|}\]

La valeur maximale est atteinte en paroi interne (\(r = R_i\)) :

\[\sigma_{VM}^{\max} = \frac{3}{2}\,\frac{p_i}{1 - (R_i/R_e)^3}\]

12.5.2 Coefficient de sécurité

\[\boxed{s = \frac{\sigma_{\max}}{\sigma_{VM}^{\max}}}\]

12.6 Renfort composite

On considère un renfort externe d’épaisseur \(t = 1\) cm en matériau composite orthotrope à fibres de carbone enroulées, avec un module tangentiel \(E_t = 250\) GPa et un module radial \(E_r = 250\) GPa. Le renfort modifie la rigidité radiale de la cuve : - les déplacements sont réduits par rapport au cas non renforcé ; - la contrainte dans la partie métallique diminue ; - la composante radiale \(\sigma_{rr}\) est continue à l’interface métal/composite ; - la composante tangentielle \(\sigma_{\theta\theta}\) présente un saut (discontinuité du matériau).

Le renfort permet d’augmenter la pression de service ou de réduire l’épaisseur de la paroi métallique pour un même niveau de sécurité.

import numpy as np
import matplotlib.pyplot as plt
from cuve_params import CuveParams

p = CuveParams()
r = np.linspace(p.R_i, p.R_e, 300)

bg = (40/255, 43/255, 48/255, 1.0)

# --- Température ---
T_r = CuveParams.T_analytique(r)

fig, ax = plt.subplots(figsize=(8, 5))
fig.patch.set_facecolor(bg); ax.set_facecolor(bg)
ax.plot((r - p.R_i)*1e3, T_r, 'tab:orange', linewidth=2)
ax.axhline(p.T_int, color='dodgerblue', ls='--', lw=0.8, alpha=0.6, label=f'$T_i = {p.T_int:.0f}$°C')
ax.axhline(p.T_ext, color='tab:red', ls='--', lw=0.8, alpha=0.6, label=f'$T_e = {p.T_ext:.0f}$°C')
ax.set_xlabel('$r - R_i$ (mm)', fontsize=14, color='grey')
ax.set_ylabel('$T$ (°C)', fontsize=14, color='grey')
ax.set_title('Température dans la paroi', fontsize=14, color='grey')
ax.legend(fontsize=11); ax.grid(True, alpha=0.3)
ax.tick_params(colors='grey')
for sp in ax.spines.values(): sp.set_color('grey')
plt.tight_layout(); plt.show()

import numpy as np
import matplotlib.pyplot as plt
from cuve_params import CuveParams

p = CuveParams()
r = np.linspace(p.R_i, p.R_e, 300)

bg = (40/255, 43/255, 48/255, 1.0)

# --- Pression interne ---
srr_p = CuveParams.sigma_rr_pression(r)
stt_p = CuveParams.sigma_tt_pression(r)
ur_p = CuveParams.u_r_pression(r)
svm_p = CuveParams.sigma_vm(srr_p, stt_p)

fig, axes = plt.subplots(1, 3, figsize=(16, 5))
fig.patch.set_facecolor(bg)

labels = [
    (r'$\sigma_{rr}$', srr_p, 'tab:green'),
    (r'$\sigma_{\theta\theta}$', stt_p, 'tab:orange'),
    (r'$\sigma_{VM}$', svm_p, 'tab:red'),
]
for ax, (label, val, color) in zip(axes, labels):
    ax.set_facecolor(bg)
    ax.plot((r - p.R_i)*1e3, val/1e6, color=color, linewidth=2)
    ax.set_xlabel('$r - R_i$ (mm)', fontsize=13, color='grey')
    ax.set_ylabel('Contrainte (MPa)', fontsize=13, color='grey')
    ax.set_title(label, fontsize=14, color='grey')
    ax.grid(True, alpha=0.3); ax.tick_params(colors='grey')
    for sp in ax.spines.values(): sp.set_color('grey')

plt.suptitle('Pression interne $p_i = 350$ bar seuls', fontsize=15, color='grey', y=1.02)
plt.tight_layout(); plt.show()

print(f"sigma_VM max = {svm_p.max()/1e6:.2f} MPa")
print(f"Coef. securite = {p.sigma_max/svm_p.max():.2f}")
print(f"u_r(R_i) = {ur_p[0]*1e3:.4f} mm")
print(f"u_r(R_e) = {ur_p[-1]*1e3:.4f} mm")

sigma_VM max = 211.11 MPa
Coef. securite = 0.95
u_r(R_i) = 0.6689 mm
u_r(R_e) = 0.6106 mm

import numpy as np
import matplotlib.pyplot as plt
from cuve_params import CuveParams

p = CuveParams()
r = np.linspace(p.R_i, p.R_e, 300)

bg = (40/255, 43/255, 48/255, 1.0)

_, axes = plt.subplots(1, 2, figsize=(14, 5))
fig = axes[0].get_figure()
fig.patch.set_facecolor(bg)

for ax in axes:
    ax.set_facecolor(bg)

# Déplacement radial
ur_p = CuveParams.u_r_pression(r)
axes[0].plot((r - p.R_i)*1e3, ur_p*1e3, 'tab:blue', linewidth=2)
axes[0].set_xlabel('$r - R_i$ (mm)', fontsize=14, color='grey')
axes[0].set_ylabel('$u_r$ (mm)', fontsize=14, color='grey')
axes[0].set_title('Déplacement radial $u_r$', fontsize=14, color='grey')
axes[0].grid(True, alpha=0.3); axes[0].tick_params(colors='grey')
for sp in axes[0].spines.values(): sp.set_color('grey')

# Contraintes
srr_p = CuveParams.sigma_rr_pression(r)
stt_p = CuveParams.sigma_tt_pression(r)
axes[1].plot((r - p.R_i)*1e3, srr_p/1e6, 'tab:green', linewidth=2, label=r'$\sigma_{rr}$')
axes[1].plot((r - p.R_i)*1e3, stt_p/1e6, 'tab:orange', linewidth=2, label=r'$\sigma_{\theta\theta}$')
axes[1].set_xlabel('$r - R_i$ (mm)', fontsize=14, color='grey')
axes[1].set_ylabel('Contrainte (MPa)', fontsize=14, color='grey')
axes[1].set_title('Contraintes dans la paroi', fontsize=14, color='grey')
axes[1].legend(fontsize=11); axes[1].grid(True, alpha=0.3); axes[1].tick_params(colors='grey')
for sp in axes[1].spines.values(): sp.set_color('grey')

plt.tight_layout(); plt.show()

12.7 Résumé des grandeurs clés

Grandeur	Formule	Valeur numérique
\(\sigma_{rr}(R_i)\)	\(-p_i\)	\(-35.00\) MPa
\(\sigma_{\theta\theta}(R_i)\)	\(p_i\frac{1 + (R_e/R_i)^3/2}{(R_e/R_i)^3 - 1}\)	\(176.11\) MPa
\(\sigma_{VM}^{\max}\)	\(|\sigma_{\theta\theta}(R_i) - \sigma_{rr}(R_i)|\)	\(211.11\) MPa
Coefficient de sécurité \(s\)	\(\sigma_{\max} / \sigma_{VM}^{\max}\)	\(0.95\)
\(u_r(R_i)\)	—	\(0.67\) mm
\(u_r(R_e)\)	—	\(0.61\) mm

Remarque : Le coefficient de sécurité \(s < 1\) indique que la pression \(p_i = 350\) bar dépasse la limite élastique du matériau pour cette géométrie. En pratique, il faudrait augmenter l’épaisseur ou utiliser un renfort composite.

12.8 Conclusion

Ce problème illustre le couplage thermo-mécanique dans une structure à symétrie sphérique : - le champ de température est solution du Laplacien sphérique (décroissance en \(1/r\)) ; - un gradient thermique induit des contraintes thermoélastiques même en l’absence de chargement mécanique extérieur ; - la solution de Lamé pour la pression interne fournit un cas-test analytique pour valider le calcul éléments finis axisymétrique ; - l’ajout d’un renfort composite modifie la rigidité et réduit les contraintes dans la partie métallique.