Données : conditions initiales {$\epsilon_n, \epsilon_n^p, (V_k)_n$} et incrément de déformation $\Delta\epsilon$

for G = 1 to N_b points de Gauss de l'élément do

Prédiction élastique : $\sigma_{n+1}^* = \sigma_n + \lambda(\text{tr } \Delta\epsilon)I + 2\mu\Delta\epsilon$

Calcul du critère : $f_{n+1}^* = \sigma_{eq_{n+1}}^* - hp_n - \sigma_y$

if $f_{n+1}^* \le 0$ then

Évolution élastique : actualiser $(\bullet)_{n+1} = (\bullet)_{n+1}^*$ & Sortie

else

Évolution élastoplastique : return mapping

while $|f_{n+1}^{(i)}| \ge \text{TOL} \quad |f_{n+1}^{(0)}| \text{ do}$

Calculer :

$$ \delta p = \frac{f_{n+1}^{(i)}}{\left. \frac{\partial f}{\partial \sigma} \right|_{n+1}^{(i)} : \mathbf{K} : \left. \frac{\partial f}{\partial \sigma} \right|_{n+1}^{(i)} - h \left. \frac{\partial f}{\partial R} \right|_{n+1}^{(i)}} $$

Corriger les variables :

$$ \sigma_{n+1}^{(i+1)} = \sigma_{n+1}^{(i)} - \delta p \mathbf{K} : \left. \frac{\partial f}{\partial \sigma} \right|_{n+1}^{(i)} $$

$$ R_{n+1}^{(i+1)} = R_{n+1}^{(i)} + h\delta p $$

$$ p_{n+1}^{(i+1)} = p_{n+1}^{(i)} + \delta p $$

$$ i \leftarrow i + 1 $$

end while

Actualiser les variables :

$$ \sigma_{n+1} = \sigma_{n+1}^{(i)} $$

$$ R_{n+1} = R_{n+1}^{(i)} $$

$$ p_{n+1} = p_{n+1}^{(i)} $$

$$ \epsilon_{n+1}^p = \epsilon_{n+1} - \mathbf{K}^{-1} : \sigma_{n+1} $$

end if

end for

**Algorithmme 4.3 – Algorithme du Return Mapping généralisé**

Dans cette méthode, aucune hypothèse n'est formulée sur l'évolution des paramètres du modèle de comportement sur le pas de temps. Les hypothèses de linéarité sont faites au niveau des corrections successives. On permet donc une discrétisation temporelle moins fine qu'avec le retour radial, tout en conservant la même précision de calcul. Aussi, au cours de la résolution, les états intermédiaires n'ont pas de sens physique mais la normalité de l'écoulement plastique à $t_{n+1}$ est assurée. Cette méthode est très attractive pour les calculs de grande taille. Cependant, un opérateur tangent ne peut être construit puisqu'une linéarisation exacte de l'algorithme ne peut être formulée. Cette méthode sera alors avantageusement combinée à des méthodes de type quasi-Newton pour le calcul de la matrice de rigidité.

Dans le cas très particulier présenté ici, on observe que la correction sur la déformation plastique cumulée (4.49) est identique à l'itération 0 à celle donnée par le retour radial pour ce même modèle (4.24), l'algorithme convergera en une itération.

La méthode du return mapping est une méthode itérative, ce qui diffère fondamentalement du retour radial ; cependant, et à ce titre, la vérification du critère à la fin du pas de temps n'est réalisée qu'à une tolérance près, définie par l'utilisateur.
Une interprétation géométrique de cet algorithme est présentée sur la figure 4.5 dans le cas d'un comportement parfaitement plastique. À chaque itération, les variables sont calculées et mises à jour par projection normale à des isopotentialles du critère des mêmes quantités calculées à l'itération précédente.

<figure><img src="image_71.png" /><figcaption>**Figure 4.5 – Return mapping**</figcaption></figure>
## A.1 Plasticité monodimensionnelle

### A.1.1 Modèle de Saint-Venant

On se propose d'étudier dans un premier temps le modèle rhéologique de Saint-Venant illustré en figure A.1.

<figure><img src="image_72.png" /><figcaption>Figure A.1 – Modèle de Saint-Venant</figcaption></figure>

1. Préciser les lois de comportement des éléments constituant ce modèle.

2. Tracer la courbe de charge monotone puis cyclique de ce modèle. À quel type de solide élastoplastique le modèle de Saint-Venant correspond-t-il ?

3. Effectuer un bilan d'énergie sur un cycle charge-décharge puis pour un cycle traction /compression fermé. Calculer les énergies stockée et dissipée associées.

### A.1.2 Plasticité en chargement monotone et cyclique

On considère ici un matériau élastoplastique avec écrouissage linéaire, soit isotrope de module K, soit cinématique de module C, tel qu'illustré en figure A.2.

<figure><img src="image_73.png" /><figcaption>Figure A.2 – Écrouissage cinématique linéaire</figcaption></figure>

1. Tracer la courbe de traction monotone pour chaque modèle. Calculer le module tangent associé à la phase élastoplastique de chaque modèle.

2. On note $\sigma_A$ la contrainte de traction maximale appliquée avant décharge élastique. Déterminer pour quelle contrainte le matériau replastifie en compression. Quelle différence y a-t-il entre ces deux modèles ?

3. Donner l'allure de la réponse de chaque modèle pour un chargement cyclique symétrique. Quel modèle vous semble le plus adapté en cyclique ?
4. Effectuer un bilan d'énergie sur un cycle charge-décharge puis pour un cycle traction /compression fermé. Calculer les énergies stockée et dissipée associées.

### A.1.3 Modèle de Saint-Venant généralisé

Le modèle de Saint-Venant généralisé, représenté sur la figure A.3, permet d'intégrer une modélisation de l'écrouissage par morceaux.

<figure><img src="image_74.png" /><figcaption>Figure A.3 – Modèle de Saint-Venant généralisé</figcaption></figure>

1. Tracer la courbe de charge monotone puis cyclique de ce modèle.

2. Effectuer un bilan d'énergie sur un cycle charge-décharge puis pour un cycle traction /compression fermé. Calculer les énergies stockée et dissipée associées.

## A.2 Loi de comportement élastoplastique

On propose de développer les équations d'un modèle élasto-plastique associé dans le cadre de la thermodynamique des processus irréversibles. L'énergie libre $\psi$ du modèle s'écrit :

$$ \rho\psi(\boldsymbol{\varepsilon}^e, \boldsymbol{\alpha}, r) = \rho\psi_e(\boldsymbol{\varepsilon}^e) + \rho\psi_{\boldsymbol{\alpha}}(\boldsymbol{\alpha}) + \rho\psi_r(r) \quad (A.1) $$

avec

$$ \rho\psi_e(\boldsymbol{\varepsilon}^e) = \frac{1}{2}\mathbf{K} : \boldsymbol{\varepsilon}^e : \boldsymbol{\varepsilon}^e \quad ; \quad \rho\psi_{\boldsymbol{\alpha}}(\boldsymbol{\alpha}) = \frac{1}{3}H\boldsymbol{\alpha} : \boldsymbol{\alpha} \quad ; \quad \rho\psi_r(r) = \frac{1}{2}Qr^2 $$

où $\mathbf{K}$ est le tenseur d'élasticité, $H$ et $Q$ sont des constantes matériau. Les variables internes associées aux écrouissages cinématique et isotrope sont respectivement notées $\boldsymbol{\alpha}$ et $r$. Les variables associées du modèle sont les couples suivants $(\boldsymbol{\sigma}, \boldsymbol{\varepsilon}^e)$, $(\mathbf{X}, \boldsymbol{\alpha})$ et $(R, r)$. On donne d'autre part l'expression du convexe d'élasticité basé sur le critère de von Mises :

$$ f(\boldsymbol{\sigma}, \boldsymbol{\alpha}, r) = \sqrt{\frac{3}{2}(\mathbf{s} - \mathbf{X}) : (\mathbf{s} - \mathbf{X}) - R - \sigma_y} \le 0 \quad (A.2) $$

où $\mathbf{s}$ est le déviateur des contraintes et $\mathbf{X}$ le tenseur (déviateur) repérant le centre du convexe d'élasticité dans l'espace des contraintes.

1. Écrire les lois d'état du modèle.

2. Écrire les lois d'évolution du modèle.

3. Calculer la vitesse de déformation plastique cumulée. Montrer que $\dot{\mathbf{p}} = \dot{\boldsymbol{\lambda}} = \dot{r}$ pour ce modèle.

4. Écrire la dissipation intrinsèque $\mathcal{D}_{\text{int}}$.

5. Écrire l'équation de cohérence, et en déduire l'expression du multiplicateur plastique.

6. Calculer l'opérateur tangent $\mathbf{H} = \frac{\partial \dot{\boldsymbol{\sigma}}}{\partial \boldsymbol{\varepsilon}}$.
7. On considère à présent un essai de traction :

(a) Tracer la courbe ($\varepsilon^p, \sigma$) pour un cycle de charge-décharge.

(b) Donner l'allure de la courbe contrainte-déformation pour le même cycle. Déterminer le module tangent.

## A.3 Torsion élastoplastique d'un arbre cylindrique plein

On considère un cylindre circulaire de rayon $R$, d'axe $(O, \mathbf{z})$ et de hauteur $h$. L'origine $O$ est placée sur la face inférieure. On note $(r, \theta, z)$ les coordonnées cylindriques et $(e_r, e_\theta, e_z)$ la base orthonormée naturellement associée.

La face inférieure est fixe et on impose à la face supérieure une rotation $\alpha h$ autour de l'axe $(O, \mathbf{z})$ où $\alpha$ est l'angle unitaire de torsion, de façon croissante au cours du temps $\dot{\alpha} > 0$. Le bord latéral est non chargé et on néglige les effets de la pesanteur. On note $C$ le couple de torsion exercé sur la face supérieure.

Le matériau est supposé élastique-parfaitement plastique et obéit au critère de von Mises et à la loi d'écoulement plastique associée par normalité. On note $\mu$ le module de cisaillement (second coefficient de Lamé) et $\sigma_0$, la limite d'élasticité (constante) en traction simple.

Le problème est traité dans le cadre des petites transformations et on suppose que le cylindre est dans son état naturel, non plastifié, à l'instant initial. On rappelle la solution purement élastique du problème en coordonnées cylindriques : $u_\theta(r, z) = \alpha r z$ et $\sigma_{\theta z}(r) = \mu \alpha r$, les autres composantes de $\mathbf{u}$ et $\sigma$ étant nulles.

1. Montrer que

$$C = \int_{0}^{R} 2\pi r^2 \sigma_{\theta z}(r) dr \qquad (A.3)$$

En déduire la valeur de $C$ en fonction du paramètre de chargement $\alpha$ pour la solution purement élastique du problème.

2. Par analogie avec la solution purement élastique, on introduit l'hypothèse que $\sigma_{\theta z}$ est aussi la seule composante non nulle dans la solution élastoplastique. Montrer que $\sigma_{\theta z}$ ne dépend que de $r$. On rappelle l'expression des équations d'équilibre en coordonnées cylindriques :

$$ \begin{aligned} \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} + \frac{\sigma_{rr} - \sigma_{\theta\theta}}{r} &= 0 \\ \frac{\partial \sigma_{\theta r}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta\theta}}{\partial \theta} + \frac{\partial \sigma_{\theta z}}{\partial z} + \frac{2\sigma_{\theta r}}{r} &= 0 \\ \frac{\partial \sigma_{zr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{z\theta}}{\partial \theta} + \frac{\partial \sigma_{zz}}{\partial z} + \frac{\sigma_{zr}}{r} &= 0 \end{aligned} \qquad (A.4) $$

3. En écrivant la loi d'élasticité et la loi d'écoulement plastique, montrer que toutes les composantes de $\dot{\mathbf{e}}^e, \dot{\mathbf{e}}^p$ et $\dot{\mathbf{e}}$ sont nulles à l'exception des composantes « $\theta z$ ». On rappelle l'expression des composantes de $\mathbf{e}$ en coordonnées cylindriques :

$$ \begin{aligned} \varepsilon_{rr} &= \frac{\partial u_r}{\partial r} & \varepsilon_{\theta\theta} &= \frac{1}{r} \left( \frac{\partial u_\theta}{\partial \theta} + u_r \right) & \varepsilon_{zz} &= \frac{\partial u_z}{\partial z} \\ \varepsilon_{r\theta} &= \frac{1}{2} \left[ \frac{1}{r} \left( \frac{\partial u_r}{\partial \theta} - u_\theta \right) + \frac{\partial u_\theta}{\partial r} \right] & \varepsilon_{\theta z} &= \frac{1}{2} \left( \frac{\partial u_\theta}{\partial z} + \frac{1}{r} \frac{\partial u_z}{\partial \theta} \right) & \varepsilon_{zr} &= \frac{1}{2} \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r} \right) \end{aligned} \qquad (A.5) $$

4. On introduit maintenant l'hypothèse (naturelle compte tenu de la symétrie de révolution) que $\mathbf{u}$ ne dépend pas de $\theta$. Montrer que les composantes $u_r$ et $u_z$ sont nulles et que $u_\theta$ est de la forme (à chargement fixé) :

$$u_{\theta}(r, z) = f(z)r \qquad (A.6)$$
5. On introduit l'hypothèse que $\boldsymbol{\epsilon}$ ne dépend pas de $z$. Montrer que :

$$u_{\theta}(r, z) = \alpha r z \qquad (A.7)$$

6. En écrivant la loi d'élasticité et la loi d'écoulement plastique pour les composantes « $\theta z$ » des déformations, et en distinguant le cas de la zone élastique et de la zone plastique, montrer que $\sigma_{\theta z} \ge 0$ partout. En déduire la valeur de $\sigma_{\theta z}$ dans la zone plastique.

7. On suppose désormais que la zone $r \le a$ est élastique et la zone $r \ge a$, plastique. Calculer $\sigma_{\theta z}$ dans la zone élastique. En se plaçant en $r = a$, en déduire l'expression de $a$ en fonction de $\alpha, \mu$ et $\sigma_0$. À quelle condition sur $\alpha$ le cylindre est-il purement élastique ?

8. On se place dans le cas élastoplastique ($a < R$). Calculer le couple de torsion $C$ en fonction de $\alpha, R \mu$ et $\sigma_0$. Indiquer les valeurs de $C^{(1)}$, la charge de première plasticité, et de $C^{(l)}$, la charge limite. Représenter graphiquement la courbe de charge $C = f(\alpha)$. À quoi l'augmentation de $C$ en fonction de $\alpha$ dans la plage comprise entre $C^{(1)}$ et $C^{(l)}$ est-elle due ?

9. On impose une torsion $\alpha_0$ suffisamment grande pour quasiment entièrement plastifier le cylindre ($a \simeq 0$), le couple $C$ atteignant donc quasiment la valeur limite $C^{(l)}$. On ramène ensuite ce couple à zéro. Calculer les contraintes résiduelles après décharge. On fera l'hypothèse, qui sera vérifiée *a posteriori*, que la décharge est purement élastique.

## A.4 Implémentation de lois de comportement à seuil

### A.4.1 Intégration numérique temporelle

1. Rappeler l'utilisation des schémas des trapèzes généralisés et des différences centrales généralisées pour la résolution discrète d'une équation différentielle ordinaire d'ordre un :

$$\dot{y} = G(y, t), \quad y(t = 0) = 0 \qquad (A.8)$$

où $G(y, t)$ est une fonction connue.

2. Spécifier les schémas précédents pour retrouver les schémas d'Euler explicite et implicite. Préciser les conditions de stabilité des algorithmes.

### A.4.2 Implémentation numérique d'un modèle de plasticité

On s'intéresse au modèle élastoplastique avec l'écrouissage cinématique non linéaire d'Armstrong-Frédéricks dont on rappelle les équations constitutives :

$$ 
\begin{aligned}
\dot{\boldsymbol{e}} &= \dot{\boldsymbol{e}}^e + \dot{\boldsymbol{e}}^p \quad ; \quad \dot{\boldsymbol{e}}^e = \boldsymbol{K}^{-1} : \dot{\boldsymbol{\sigma}} \quad ; \quad \dot{\boldsymbol{e}}^p = \frac{3}{2} \dot{\boldsymbol{\lambda}} \frac{\boldsymbol{s} - \boldsymbol{X}}{(\boldsymbol{s} - \boldsymbol{X})_{\text{eq}}} \\
\dot{\boldsymbol{X}} &= \frac{2}{3} C \dot{\boldsymbol{e}}^p - \gamma \boldsymbol{X} \dot{\boldsymbol{p}} \\
\dot{\boldsymbol{\lambda}} &= \dot{\boldsymbol{p}} = \sqrt{\frac{2}{3} \dot{\boldsymbol{e}}^p : \dot{\boldsymbol{e}}^p} \quad ; \quad f = (\boldsymbol{s} - \boldsymbol{X})_{\text{eq}} - \sigma_y \le 0 \\
\dot{\boldsymbol{\lambda}} \ge 0, \quad f \le 0, \quad \dot{\boldsymbol{\lambda}} f = 0
\end{aligned}
\qquad (A.9) $$

Ce système comprend la partition du taux de déformation, la loi d'élasticité, la loi d'écoulement, la loi d'évolution de l'écrouissage cinématique, l'expression de la déformation plastique cumulée, le critère basé sur la contrainte équivalente de von Mises et les conditions de complémentarité de Kühn-Tucker.

1. Appliquer le schéma d'Euler implicite aux équations du modèle de plasticité considéré.
<figure><img src="image_75.png" /></figure>

2. Écrire la prédiction élastique à l'instant $t_{n+1}$, en supposant connu d'une part l'état mécanique à l'instant $t_n$ ($\epsilon_n, \epsilon_n^p, \sigma_n, X_n$), et d'autre part la déformation totale $\epsilon_{n+1} = \epsilon_n + \Delta\epsilon$ à l'instant $t_{n+1}$.

3. Écrire le test à effectuer sur le critère $f$ et ses conséquences quant à la résolution numérique.
Poser le problème comme un système d'équation non linéaire à résoudre.

4. Proposer au moins deux méthodes numériques pour résoudre ce système d'équation non linéaire. Discuter.

5. Appliquer la méthode de Newton au modèle considéré. Calculer le module tangent cohérent.

6. Préciser les tests de convergence possibles.
## B.1 Sphère sous pression en élastoplasticité [2012-2013]

On considère un récipient sphérique creux (figure C.2) de rayon intérieur $a$, de rayon extérieur $b$, soumis à une pression interne $p$ croissante à partir de zéro. On note $(r, \theta, \varphi)$ les coordonnées sphériques et $(e_r, e_\theta, e_\varphi)$ la base orthonormée naturellement associée. Le matériau considéré est élastique parfaitement plastique et obéit au critère de von Misès. On note $E$ le module d'Young, $\nu$ le coefficient de Poisson et $\sigma_0$ la limite d'élasticité en traction simple. Le problème est ici traité dans le cadre des petites transformations isothermes et quasi-statiques.

<figure><img src="image_76.png" /><figcaption>Figure B.1 – Sphère sous pression</figcaption></figure>

On rappelle la solution purement élastique de ce problème :

$$ 
\begin{align} 
\sigma_{rr} &= p \frac{a^3}{b^3 - a^3} \left(1 - \frac{b^3}{r^3}\right) \nonumber \\ 
\sigma_{\theta\theta} &= \sigma_{\varphi\varphi} = p \frac{a^3}{b^3 - a^3} \left(1 + \frac{b^3}{2r^3}\right) \tag{B.1} \\ 
u_r &= \frac{p}{E} \frac{a^3}{b^3 - a^3} \left((1 - 2\nu)r + (1 + \nu)\frac{b^3}{2r^2}\right) \nonumber 
\end{align} 
$$

On admettra que dans le cas élastoplastique, les seules composantes non nulles du tenseur des contraintes sont encore $\sigma_{rr}$, $\sigma_{\theta\theta}$ et $\sigma_{\varphi\varphi}$, et que ces composantes ne dépendent que du rayon.

1. **Fin de la phase élastique.** À partir de l'expression de la solution purement élastique (B.1), calculer la contrainte équivalente de von Misès (on utilisera avantageusement le fait que le critère de Misès est indépendant de la pression hydrostatique). En appliquant le critère du même nom, en déduire la charge $p^{(1)}$ de première plastification de la structure. Où la plasticité apparaît-elle en premier dans la structure ?
2. **Phase élastoplastique.** On suppose à présent que $p > p^{(1)}$, et que la structure est composée d'une zone élastique et d'une zone plastique, délimitées par le rayon $c$ tel que $a \le c \le b$, dont on précisera les domaines respectifs.

(a) Donner la solution en contraintes dans la zone élastique en fonction de l'inconnue auxiliaire $p' = -\sigma_{rr}(r = c)$.

(b) En écrivant que le critère est atteint en $r = c$, déterminer l'inconnue auxiliaire $p'$, et ré-exprimer la solution en contraintes dans la zone élastique en fonction du rayon plastique $c$.

(c) On considère maintenant la zone plastique. On rappelle que l'équation d'équilibre radiale s'écrit :

$$ \frac{d\sigma_{rr}}{dr} + \frac{2(\sigma_{rr} - \sigma_{\theta\theta})}{r} = 0 \qquad (B.2) $$

En combinant l'équation d'équilibre radial (B.2) et le critère (supposé atteint dans la zone plastique) en admettant que le signe de $\sigma_{rr} - \sigma_{\theta\theta}$ est le même dans la zone considérée, déterminer la forme analytique de $\sigma_{rr}$.

(d) En écrivant une condition de raccord convenable (que l'on précisera) des solutions des zones élastique et plastique, calculer le champ de contraintes dans la zone plastique. Tracer l'évolution des composantes de contraintes en fonction du rayon de la sphère. Commenter la continuité des composantes de contraintes et l'équilibre (B.2) le long du rayon. Pourquoi les composantes $\sigma_{\theta\theta} = \sigma_{\varphi\varphi}$ sont-elles continues en $r = c$ ?

(e) En écrivant la condition aux limites en $r = a$, calculer $p$ en fonction de $c$. Montrer que l'expression obtenue définit $c$ comme une fonction croissante de $p$.

(f) Déterminer la valeur de la charge limite $p^{(L)}$. Vous préciserez ce qu'est une charge limite et la raison pour laquelle il existe ou pas une charge limite pour une structure élastoplastique. Tracer la courbe de charge de la structure. Pourquoi la pression augmente t-elle dans la phase élastoplastique de la structure ? À quoi cela est-il dû ?

3. **Décharge à partir de la charge limite.** On suppose dans cette partie qu'on a chargé la structure jusqu'à la charge limite, puis qu'on décharge complètement jusqu'à $p = 0$. Afin de regrouper les éventuels décrocheurs, on donne l'expression des contraintes dans la sphère à la charge limite :

$$ \sigma_{rr} = 2\sigma_0 \ln\left(\frac{r}{b}\right) \qquad (B.3) \\ \sigma_{\theta\theta} = \sigma_{\varphi\varphi} = \sigma_{rr} + \sigma_0 $$

(a) En supposant la décharge purement élastique, déterminer les variations $\Delta\sigma_{rr}$, $\Delta\sigma_{\theta\theta}$ et $\Delta\sigma_{\varphi\varphi}$ des contraintes pendant la décharge, puis la valeur de ces contraintes après décharge.

(b) Déterminer en quel point peut éventuellement commencer la recharge plastique pendant la décharge, et à quelle condition sur $a$ et $b$ cette éventualité ne se produit pas (c'est-à-dire que la décharge est effectivement purement élastique). On posera le paramètre $\lambda = \frac{a^3}{b^3}$, et on pourra exprimer la condition de décharge élastique sous la forme $\lambda > \lambda_0$, $\lambda_0$ étant à déterminer.

## B.2 Analyse d'un système de barres en traction en élastoplasticité [2013-2014]

On considère le système de barres en traction de la figure B.2 constituée de trois barres de longueurs identiques $L$ et de sections identiques $S$, articulées par rapport au bâti et par rapport à
<figure><img src="image_77.png" /></figure>

la poutre ABC supposée indéformable. Les distances AB et BC sont aussi supposées identiques et égales à L.

Ce système est soumis à une force verticale descendante P concentrée au point A. On suppose que les barres sont à l'état initial sans contraintes, et que le chargement est appliqué suffisamment lentement. On suppose d'autre part que les trois barres ont un comportement élastique parfaitement plastique identique et obéissent au critère de von Misès. On note E, ν et σ₀ le module d'Young, le coefficient de Poisson et la limite d'élasticité en traction des barres respectivement. Le problème est ici traité dans le cadre des petites transformations.

<figure><img src="image_78.png" /><figcaption>Figure B.2 – Système de trois barres en traction</figcaption></figure>

On se propose ici d'étudier cette structure numériquement avec le code de calcul CAST3M, et de comparer la solution numérique avec la solution analytique existant pour ce problème.

## Partie 1 : structure purement élastique

1. Faire le maillage des trois barres 1, 2 et 3, en n'utilisant qu'un seul élément fini "barre" par barre. La poutre indéformable ne sera pas maillée, son influence sera prise en compte au travers d'une condition aux limites. On prendra $L = 100$ mm, une charge P de 20 kN, des barres en acier de construction classique S200 ($E = 210\,000$ MPa, $\nu = 0.3$, $\sigma_0 = 200$ MPa) de diamètre 10 mm.

2. En isolant la poutre indéformable, montrer que le système décrit sur la figure B.2 est hyperstatique de degré 1.

3. La présence de la poutre indéformable ABC implique une condition de compatibilité géométrique entre les déplacements verticaux des points A, B et C ; en notant $\delta_i$ l'allongement de la $i$-e barre, et en raisonnant sur la configuration déformée du système, montrer que cette condition s'écrit :

$$ \delta_1 + \delta_3 - 2\delta_2 = 0 \qquad (B.4) $$

4. L'opérateur RELA dans CAST3M permet d'imposer des relations entre des degrés de liberté définis aux noeuds du maillage. Dans notre cas, la syntaxe prend la forme suivante :

RIG1 = RELA VAL1 MOT1 GEO1 + VAL2 MOT2 GEO2 + ...

où VAL1 est un coefficient multiplicatif de MOT1 qui décrit la composante d'intérêt ('UY' ici) du déplacement de l'entité géométrique décrite par GEO1.

Résolvez de façon numérique le problème purement élastique, et comparer avec la solution analytique :

$$ N_1 = \frac{5P}{6}; N_2 = \frac{P}{3}; N_3 = -\frac{P}{6}; \delta_i = \frac{N_i L}{ES}, \quad i = 1, 2, 3 \qquad (B.5) $$
où les $N_i$ ($i = 1, 2, 3$) sont les efforts normaux dans les barres.

## Partie 2 : structure élastique parfaitement plastique

5. Quelle barre plastifiera en premier ? À partir de la solution purement élastique (B.5), donner la charge $P^{(1)}$ de première plastification de la structure.

6. Définissez un trajet de chargement monotone croissant dans CAST3M, et définissez les propriétés nécessaires dans une table pour l'opérateur PASAPAS.

7. On note $\Delta P$ l'incrément de charge entre la charge de première plastification $P^{(1)}$ et la charge appliquée. La solution analytique élastoplastique donne l'allongement de la barre 1 :

$$ \delta_1 = \delta_1^{(1)} + \frac{L}{ES} 5\Delta P \qquad (B.6) $$

où $\delta_1^{(1)}$ est l'allongement de la barre 1 à la charge de première plastification $P^{(1)}$. Comparer les solutions analytique et numérique. Commenter.

8. Que se passe t-il si on augmente la charge $P$ de 20 à 25 kN, et pourquoi ? À quoi ce phénomène correspond-il ?

9. On donne le chargement $P = 23 550$ N. Faites une décharge à partir de l'état chargé avec cette dernière charge. La décharge est-elle purement élastique ? Comparer l'état de contrainte résiduel dans les barres avec la solution analytique :

$$ N_1^{\text{res}} = -\frac{N_0}{4}; \quad N_2^{\text{res}} = \frac{N_0}{2}; \quad N_3^{\text{res}} = -\frac{N_0}{4} \qquad (B.7) $$

où $N_0 = \sigma_0 S$ est l'effort associé à la limite d'élasticité en traction/compression des barres. Commenter.

## B.3 Implémentation d'une loi de comportement élastoplastique [2014-2015]

Le rendu de votre devoir se composera d'une copie, sur laquelle vous rédigerez vos réponses aux questions, et de fichiers numériques que vous posterez directement sur le serveur pédagogique dans le travail "Examen DEFIR 2015", à votre nom. L'évaluation est individuelle. Les fichiers MatLab ainsi que les figures (au format .pdf, .eps, .ps ou .png) associées aux tracés demandés seront téléchargés sur le serveur avant la fin de l'examen. Attention, rien ne peut plus être téléchargé après 12h15.

**Commentaire** Seuls les documents accessibles sur le serveur pédagogiques sont autorisés pendant l'examen. Tout autre document n'est pas autorisé et doit être rangé.

On s'intéresse au modèle élastoplastique à écrouissage cinématique non linéaire d'Armstrong-Fredericks en petites transformations, dont on rappelle les équations constitutives en 1D :

$$ \dot{\epsilon} = \dot{\epsilon}^e + \dot{\epsilon}^p \qquad (B.8) $$

$$ \dot{\sigma} = E \dot{\epsilon}^e \qquad (B.9) $$

$$ \dot{\epsilon}^p = \dot{p} \operatorname{sgn}(\sigma - X) \qquad (B.10) $$

$$ \dot{X} = C \dot{\epsilon}^p - \gamma X \dot{p} \qquad (B.11) $$

$$ f = |\sigma - X| - \sigma_y \le 0 \qquad (B.12) $$

$$ \dot{p} \ge 0; \quad f \le 0; \quad f \dot{p} = 0 \qquad (B.13) $$
<figure><img src="image_79.png" /></figure>

où les équations (B.8), (B.9), (B.10), (B.11), (B.12), et (B.13) réfèrent à la partition du taux de déformation, la loi d'élasticité, la loi d'écoulement plastique, la loi d'évolution de la variable associée à l'écrouissage cinématique, au critère et aux conditions de complémentarité de Kuhn-Tucker respectivement. Les constantes E, C et γ sont des paramètres matériaux qui dénotent le module d'Young et les paramètres de l'écrouissage cinématique respectivement.

On se propose ici de discrétiser ce modèle de comportement par un schéma d'Euler implicite, et de l'implémenter dans un code de calcul par éléments finis.

1. Discrétiser le système d'équations (B.8)-(B.13) en appliquant le schéma d'Euler implique.

2. Un schéma de prédiction/correction est adopté pour l'intégration temporelle de cette loi de comportement. On suppose dans un premier temps un état auxiliaire : l'état élastique prédit. En notant (●)* les quantités élastiquement prédites, écrire les équations associées à cet état élastique prédit.

3. Si le critère est violé, il est nécessaire d'effectuer une correction plastique. Montrer que le système discret d'équations à résoudre s'écrit :

$$ \sigma_{n+1} = \sigma^* - E \Delta p \operatorname{sgn}(\sigma_{n+1} - X_{n+1}) \quad (B.14) $$

$$ \varepsilon_{n+1}^p = \varepsilon_n^p + \Delta p \operatorname{sgn}(\sigma_{n+1} - X_{n+1}) \quad (B.15) $$

$$ X_{n+1}(1 + \gamma \Delta p) = X_n + C \Delta p \operatorname{sgn}(\sigma_{n+1} - X_{n+1}) \quad (B.16) $$

$$ |\sigma_{n+1} - X_{n+1}| - \sigma_y = 0 \quad (B.17) $$

4. On peut observer qu'un découplage peut être opéré dans le système (B.14)-(B.17), entre les équations (B.14), (B.16), (B.17) d'une part, et (B.15) d'autre part. En effet, $\varepsilon_{n+1}^p$ n'intervient pas dans le premier groupe d'équations, elle peut être mise à jour a posteriori. En utilisant MatLab, faites un programme intégrant cette loi de comportement pour un point matériel :

$$ [S,X,EP,P] = \text{armstrongFredericks}(\text{eps}, \text{Sn}, \text{Xn}, \text{EPn}, \text{Pn}, \text{E}, \text{Sy}, \text{C}, g) $$

prenant en argument la déformation totale (qu'on pilotera), l'état mécanique $(\sigma_n, X_n, \varepsilon_n^p, p_n)$ à l'instant précédent $t_n$, et les paramètres matériaux $(E, C, \gamma)$. Ce programme retourne l'état mécanique actualisé $(\sigma_{n+1}, X_{n+1}, \varepsilon_{n+1}^p, p_{n+1})$.

Pour la correction plastique, le système d'équations (B.14), (B.16), (B.17) sera résolu à l'aide de la fonction `fsolve` sur le vecteur inconnu $\mathbf{x}_{n+1}^T = \{\sigma, X, p\}_{n+1}$ en utilisant l'algorithme de Levenberg-Marquardt, selon la syntaxe suivante :

```
options = optimset('Display','off','TolFun',1e-20, ... 'Algorithm','levenberg-marquardt') ;
x = fsolve(@x) [eq1(x(1),x(2),x(3)),eq2(x(1),x(2),x(3))],
... eq3(x(1),x(2),x(3))],x0, options) ;
```

où $\sigma \equiv x(1), X \equiv x(2), p \equiv x(3), \text{eq1}(\dots)$ réfère à la première équation à coder.

5. Tracer la courbe contrainte-déformation en utilisant les valeurs résumées dans le tableau B.1 ci-dessous.

6. Montrer que le module tangent discret défini par :

$$ H_{n+1} = \frac{\partial \sigma_{n+1}}{\partial \Delta \varepsilon} \quad (B.18) $$

s'écrit :

$$ H_{n+1} = E \left( 1 - \frac{E(1 + \gamma \Delta p)^2 \operatorname{sgn}(\sigma_{n+1} - X_{n+1})}{C + E(1 + \gamma \Delta p)^2 - \frac{\gamma X_n}{\operatorname{sgn}(\sigma_{n+1} - X_{n+1})}} \right) \quad (B.19) $$
<table><caption>**Tableau B.1 – Données numériques**</caption><thead><tr><td>Paramètre</td><td>Valeur</td></tr></thead><tbody><tr><td>ε</td><td>∈ [0, 0.5]</td></tr><tr><td>E</td><td>110 GPa</td></tr><tr><td>σ<sub>y</sub></td><td>955 MPa</td></tr><tr><td>C</td><td>1000 MPa</td></tr><tr><td>γ</td><td>1</td></tr></tbody></table>

**Commentaire** Vous ne passerez pas plus de 10 minutes sur cette question.

7. On considère à présent la pâle de turbine en régime stationnaire que vous avez étudié en projet numérique, modélisée par un système unidimensionnel et soumise à un chargement centrifuge (voir figure C.1). On souhaite à présent intégrer cette loi de comportement

<figure><img src="image_80.png" /><figcaption>**Figure B.3 – Modèle unidimensionnel d'une pâle soumise à un chargement centrifuge**</figcaption></figure>

à un code de calcul par éléments finis 1D. À cet effet, le programme `bar.m` vous est fourni, vous pouvez le télécharger sur le serveur dans

Documents/ExamenDEFIR2105/programmes matlab/

Ce programme résout numériquement par la méthodes des éléments finis ce problème. La loi de comportement est appelée par la routine :

[sig_k, X_k, epsp_k, p_k, tangent_modulus] = armstrong_fredericks(x', u(:,i-1), u_k, sigma(:,i-1), X(:,i-1), epsp(:,i-1), p(:,i-1), E, Sy, C, g);

Modifiez le programme précédent que vous venez d'écrire en y intégrant une boucle sur tous les points d'intégration du maillage, et codez l'expression du module tangent (B.19). Lancez le calcul éléments finis, et tracez le déplacement, la déformation plastique ainsi que la contrainte le long de la barre pour les différents pas de chargement codés dans le programme.
**Compte rendu** Un compte rendu sera rendu à la fin du projet pour chacun des binômes. Celui-ci sera rédigé numériquement, édité au format pdf, et devra contenir les réponses aux questions rédigées dans ce texte, les tracés des courbes demandées, et la discussion des résultats et leur interprétation. Par ailleurs, vos programmes MatLab viendront compléter le dossier rendu.

## Introduction

L'objet de ce projet numérique est d'implémenter des lois de comportement élastoplastiques dans un petit code de calcul de structure par éléments finis dédié à l'analyse de structures unidimensionnelle et bidimensionnelle axisymétrique.

Deux problèmes d'évolution respectivement unidimensionnel et bidimensionnel axisymétrique seront considérés dans ce projet. Ces deux problèmes serviront de support pour comparer les résultats obtenus avec votre code aux solutions obtenues numériquement avec le code de calcul industriel par éléments finis CAST3M, et à la solution analytique pour le second problème d'évolution.

## C.1 Système unidimensionnel

### C.1.1 Développement du code de calcul par éléments finis

**Formulation du problème éléments finis 1D et algorithme** La résolution du problème d'évolution se fait ici dans le cadre des petites perturbations ; la géométrie considérée pour tous les pas de charge reste la configuration initiale du système. À chaque pas de temps (ou chaque incrément de charge en plasticité), le problème à résoudre s'écrit :

$$f^{\text{ext}} - f^{\text{int}} = 0 \qquad (C.1)$$

où $f^{\text{ext}}$ et $f^{\text{int}}$ sont respectivement les forces externes et internes, définies en une dimension telles que :

$$f^{\text{ext}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{N}^T f_d \, d\Omega + \int_{\partial \Omega_e \cap \partial \Omega_F} \mathbf{N}^T F_d \, dS \right) \qquad (C.2)$$

$$f^{\text{int}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{N}_{,x}^T \sigma \, d\Omega \right) \qquad (C.3)$$
Le système d'équations algébriques (C.1) est non linéaire, et sa résolution se fait de façon itérative. À chaque itération du calcul, le système linéaire suivant est résolu :

$$
\mathbf{K}^{(k)} \delta \mathbf{u}^{(k)} = \mathbf{R}^{(k)} \tag{C.4}
$$

avec

$$
\mathbf{K}^{(k)} = - \left. \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right|^{(k)} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{N}_{,x}^{\mathrm{T}} \mathbf{H}^{(k)} \mathbf{N}_{,x} \mathrm{d}\Omega \right) \quad \text{et} \quad \mathbf{R}^{(k)} = \mathbf{f}^{\mathrm{ext}} - (\mathbf{f}^{\mathrm{int}})^{(k)}
$$

où $\mathbf{K}^{(k)}$ est la matrice tangente calculée à l'itération $k$ du calcul, $\mathbf{H}^{(k)}$ l'opérateur tangent calculé lors de l'intégration de la loi de comportement élastoplastique. Ceci permet d'obtenir une nouvelle estimation du champ de déplacement solution de (C.1) :

$$
\mathbf{u}^{(k+1)} = \mathbf{u}^{(k)} + \delta \mathbf{u}^{(k)} \tag{C.5}
$$

Pour résoudre ce problème à chaque incrément de charge, votre programme doit suivre l’algo-
rithme C.1.

1: Faire le maillage éléments finis

2: Définition des paramètres matériaux

3: Initialisation des vecteurs (u, σ, ε, H, etc), définition de l'historique de chargement.

4: **for** *n* = 1 **to** *M* **pas de chargement do**

5: Prédiction élastique :

6: assemblage de $\mathbf{K}_{\text{elas}}$ et de $\mathbf{f}^{\text{ext}}$

7: application des conditions aux limites

8: Résolution : $\mathbf{u}^{(0)} = \mathbf{K}_{\text{elas}}^{-1} \mathbf{f}^{\text{ext}}$

9: *Intégration de la loi de comportement, pour chaque point de Gauss de chaque élément : calcul*
*de (σ, ε<sup>P</sup>, H, A<sub>k</sub>, V<sub>k</sub>)<sup>(0)</sup>*

10: Calcul des forces internes $\mathbf{f}^{\text{int}}(\mathbf{u}^{(0)})$ selon (C.3)

11: Calcul du résidu $\mathbf{R}^{(0)} = \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}(\mathbf{u}^{(0)})$ et de sa norme $\|\mathbf{R}^{(0)}\|$

12: while $||\mathbf{R}^{(k)}|| > \text{TOL do}$

13: étape globale :

14: assemblage de $\mathbf{K}^{(k)}$ selon (C.4)

15: application des conditions aux limites

16: résolution de $\mathbf{K}^{(k)} \delta \mathbf{u}^{(k)} = \mathbf{R}^{(k)}$

17: Mise à jour de $\mathbf{u}^{(k+1)}$ avec (C.5)

18: $k \leftarrow k + 1$

19: *Intégration de la loi de comportement, pour chaque point de Gauss de chaque élément :*
*calcul de (σ, ε<sup>P</sup>, H, A<sub>k</sub>, V<sub>k</sub>)<sup>(k)</sup>*

20: Calcul des forces internes $\mathbf{f}^{\text{int}}(\mathbf{u}^{(k)})$ selon (C.3)

21: Calcul du résidu $\mathbf{R}^{(k)} = \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}(\mathbf{u}^{(k)})$ et de sa norme $\|\mathbf{R}^{(k)}\|$

22: end while

23: end for

**Algorithme C.1 – Résolution du problème d’évolution élastoplastique en une dimension par la**
**méthode des éléments finis**
<figure><img src="image_81.png" /></figure>

**Intégration de la loi de comportement élastoplastique en une dimension** La loi de comportement élastoplastique particularisé en une dimension et à un écrouissage isotrope s'écrit :

$$ \epsilon = \epsilon^e + \epsilon^p \qquad (C.6) $$

$$ \sigma = E\epsilon^e \qquad (C.7) $$

$$ f = |\sigma| - R - \sigma_y \le 0 \qquad (C.8) $$

$$ \dot{\epsilon}^p = \dot{\lambda} \operatorname{sgn}(\sigma) \qquad (C.9) $$

$$ R = R(p) \qquad (C.10) $$

$$ \dot{\lambda} \ge 0; f \le 0; \dot{\lambda} f = 0 \qquad (C.11) $$

$$ \dot{\lambda} \dot{f} = 0 \qquad (C.12) $$

$$ \dot{p} = |\dot{\epsilon}^p| = \dot{\lambda} \qquad (C.13) $$

L'intégration de la loi de comportement élastoplastique est un problème local défini en chaque point d'intégration du maillage et consiste à calculer l'état mécanique {$\sigma, \epsilon^p, A_k, V_k$}$_{n+1}$ à un pas de charge, étant donné cet état au pas précédent {$\sigma, \epsilon^p, A_k, V_k$}$_n$ et un incrément de déformation totale $\Delta\epsilon$. Les algorithmes d'intégration de ce type de loi à seuil sont généralement formés d'une étape de prédiction et d'une étape de correction plastique. Les algorithmes de retour radial et de Return Mapping (algorithmes 2 et 3 du polycopié de cours) sont couramment utilisés dans les codes de calcul. en une dimension, le retour radial se simplifie comme spécifié dans l'algorithme C.2. Le module tangent $H$ se calcule en une dimension de la façon suivante :

1: Données {$\Delta\epsilon, \epsilon_n^p, (V_k)_n$}

2: for $G = 1$ to $Nb$ points de Gauss du maillage do

3: Prédiction élastique : $\sigma^* = \sigma_n + E\Delta\epsilon$

4: Calcul du critère : $f^* = |\sigma^*| - R_n - \sigma_y$

5: if $f^* \le 0$ then

6: Évolution élastique : actualiser $(\bullet)_{n+1} = (\bullet)^*$ & EXIT

7: else

8: Évolution élastoplastique : retour radial

9: Résoudre sur $\Delta p$ l'équation $\sigma^* - R(p_n + \Delta p) - \sigma_y = 0$

10: Actualiser les variables :

11: $\sigma_{n+1} = \sigma^* - E \Delta p$ signe($\sigma^*$)

12: $\epsilon_{n+1}^p = \epsilon_n^p + \Delta p$ signe($\sigma^*$)

13: $p_{n+1} = p_n + \Delta p$

14: $R_{n+1} = R(p_{n+1})$

15: Calcul du module tangent $H_{n+1}(p_{n+1})$

16: end if

17: end for

<caption>**Algorithme C.2 – Méthode du retour radial en une dimension**</caption>

$$ H = \frac{\partial \sigma}{\partial \Delta \epsilon} = E \left( 1 - \frac{\partial \Delta p}{\partial \Delta \epsilon} \operatorname{signe}(\sigma^*) \right) \qquad (C.14) $$

Il reste alors à calculer le terme $\partial\Delta p/\partial\Delta\epsilon$ en fonction du modèle d'écrouissage isotrope retenu. Celui-ci s'obtient en supposant une évolution élastoplastique $f_{n+1} = 0$ en calculant $\partial f_{n+1}/\partial\Delta\epsilon = 0$.

**Application des conditions aux limites** Plusieurs méthodes peuvent être utilisées pour appliquer les conditions aux limites essentielles dans la méthode des éléments finis. Les trois
méthodes les plus utilisées sont l'élimination, la pénalité et la méthode des multiplicateurs de Lagrange.

— l'élimination consiste à regrouper les équations du système linéaire associées aux degrés de liberté (indicées 1) et celles associées aux déplacements nodaux imposés (indicées 2). Le système linéaire se sous-divise alors de la façon suivante :

$$ \begin{pmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} \\ \mathbf{K}_{21} & \mathbf{K}_{22} \end{pmatrix} \begin{pmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \end{pmatrix} = \begin{pmatrix} \mathbf{F}_1 \\ \mathbf{F}_2 \end{pmatrix} \qquad (C.15) $$

où $\mathbf{u}_2 = \mathbf{u}^{(d)}$ sont les déplacements imposés. À chaque itération du calcul, on résoud :

$$ \mathbf{K}_{11}\mathbf{u}_1 = \mathbf{F}_1 - \mathbf{K}_{12}\mathbf{u}_2 \qquad (C.16) $$

Les équations supplémentaires sont utilisées pour déterminer les forces de réaction $\mathbf{F}_2$ si nécessaire.

— la pénalité consiste à ajouter au premier et second membre un coefficient de pénalité $k$ de valeur importante sur la ligne du degré de liberté à imposer. Si $\mathbf{1}_Q$ représente un vecteur colonne où la $Q$-ième composante est définie à 1, les autres à zéro, le système à résoudre s'écrit :

$$ (\mathbf{K} + k\mathbf{1}_Q\mathbf{1}_Q^T)\mathbf{u} = \mathbf{F} + k u_Q^{(d)}\mathbf{1}_Q \qquad (C.17) $$

Quand $k \to \infty$, $u_Q \to u_Q^{(d)}$. Attention cependant au conditionnement de la matrice qui explose pour $k$ très grand.

— la méthode des multiplicateurs de Lagrange consiste à imposer les contraintes suivantes $\mathbf{B}u = \mathbf{u}^{(d)}$; celles-ci viennent augmenter la taille du système linéaire. Les multiplicateurs de Lagrange introduit $\lambda$ représentent alors les forces nécessaires pour imposer les déplacements souhaités :

$$ \begin{pmatrix} \mathbf{K} & \mathbf{B}^T \\ \mathbf{B} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{u} \\ \mathbf{\lambda} \end{pmatrix} = \begin{pmatrix} \mathbf{F} \\ \mathbf{u}^{(d)} \end{pmatrix} \qquad (C.18) $$

il est à noter qu'il est nécessaire d'utiliser une méthode directe avec pivotage pour résoudre ce système ou une méthode de type Uzawa.

## Outils MatLab pour la rédaction du code éléments finis

★ Un maillage unidimensionnel uniforme peut être réalisé à l'aide de la fonction `linspace`, de la façon suivante :

$$ x = \text{linspace}(x_{\text{debut}}, x_{\text{fin}}, Nb_{\text{element}}) \qquad (C.19) $$

où $x$ est un vecteur contenant la position des nœuds du maillage 1D.

★ Deux routines MatLab `Matrix1D.m` et `Vector1D.m` sont mises à votre disposition pour vous faciliter le calcul des premiers et seconds membres en une dimension. Elles permettent respectivement de calculer à partir du maillage unidimensionnel $x$ une matrice associée à une forme bilinéaire $a(u^*, u)$ et un vecteur associé à une forme linéaire $l(u^*)$. L'intégration des quantités élémentaires est réalisée sur chaque élément par la méthode de Gauss, et leur assemblage se fait à la volée. Les matrice et vecteur associés à ces formes bilinéaire et linéaire font intervenir les $p$-ième et $q$-ième dérivées des fonctions de forme $N$ par rapport à la coordonnée longitudinale, ainsi que des champs soit définis aux
nœuds (nodeField) ou constants par éléments (elemField), de sorte que :

$$a(u^*, u) = \int_{\Omega} \frac{d^{(p)}u^*}{dx^{(p)}} (\text{nodeField} \times \text{elemField}) \frac{d^{(q)}u}{dx^{(q)}} d\Omega \approx u^* M u \quad (C.20)$$

$$l(u^*) = \int_{\Omega} \frac{d^{(p)}u^*}{dx^{(p)}} (\text{nodeField} \times \text{elemField}) d\Omega \approx u^* f \quad (C.21)$$

La syntaxe est alors la suivante :

$$\mathbf{M} = \text{Matrix1D}(x, p, q, \text{nodeField}, \text{elemField}) \\ = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \frac{d^{(p)}\mathbf{N}}{dx^{(p)}} (\text{nodeField} \times \text{elemField}) \frac{d^{(q)}\mathbf{N}^T}{dx^{(q)}} d\Omega \right) \quad (C.22)$$

$$\mathbf{f} = \text{Vector1D}(x, p, \text{nodeField}, \text{elemField}) \\ = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \frac{d^{(p)}\mathbf{N}^T}{dx^{(p)}} (\text{nodeField} \times \text{elemField}) d\Omega \right) \quad (C.23)$$

Pour une matrice de masse, on définira $p = q = 0$, et pour une matrice de raideur $p = q = 1$. Pour le vecteur des forces extérieures, on définira $p = 0$ et pour le vecteur des forces intérieures $p = 1$.

* Le calcul de la norme $p$ d'un vecteur $x$ se fait avec l'opérateur $\text{norm}(x, p)$.

* La résolution d'une équation ou d'un système d'équations non linéaires peut être effectuée avec la fonction `fsolve` qui résoud $g(x = 0)$.

### C.1.2 Application : pâle de turbine soumise à un chargement centrifuge

On considère une pâle de turbine dans un turboréacteur en régime stationnaire. Cette pâle est modélisée par un système unidimensionnel soumis à un chargement centrifuge modélisé par des forces linéiques variant linéairement avec le rayon courant $r$ :

$$\mathbf{f}(r) = \rho S r \omega^2 \mathbf{e}_r \quad (C.24)$$

où $\rho$, $S$ et $\omega$ sont ici la masse volumique, la section et la vitesse de rotation de la pâle respectivement. La figure C.1 présente la modélisation proposée. Cette pâle est supposée

<figure><img src="image_82.png" /><figcaption>**Figure C.1** – Modèle unidimensionnel d'une pâle soumise à un chargement centrifuge</figcaption></figure>

constituée de l'alliage de titane TA6V. On suppose que son comportement dans la gamme de vitesse considérée est élastoplastique avec un écrouissage isotrope variant en loi de puissance :

$$R(p) = Bp^n \quad (C.25)$$

où $B$ et $n$ sont des paramètres matériaux. On rappelle que $p$ est la déformation plastique cumulée. Les valeurs des paramètres matériaux de cet alliage sont résumés dans le tableau C.1. Le fichier `pale.dgibi` de la mise en données CAST3M vous est fourni et sert ici de référence pour la comparaison de vos résultats. Lors de son exécution, les données nécessaires pour la comparaison avec vos résultats sont extraits dans des fichiers textes.
<table><caption>**Tableau C.1 – Paramètres élastoplastiques du TA6V**</caption><thead><tr><td>paramètres</td><td>valeurs numériques</td></tr></thead><tbody><tr><td>ρ</td><td>4430 kg m<sup>-3</sup></td></tr><tr><td>E</td><td>110000 MPa</td></tr><tr><td>σ<sub>y</sub></td><td>955 MPa</td></tr><tr><td>B</td><td>770 MPa</td></tr><tr><td>n</td><td>0.557</td></tr></tbody></table>

1. Faire le maillage du domaine de calcul. On prendra une longueur $L = 0,2 \text{ m}$, et une section rectangulaire de $50 \times 20 \text{ mm}$ supposée constante. À la vue de l'expression du chargement centrifuge (C.24), discuter d'un maillage adapté pour une bonne représentation de la solution.

2. On considère que la pâle tourne à 1500 tr/min. Comparer dans un premier temps les solutions numériques élastiques (celle de CAST3M et la votre) à la solution analytique suivante :

$$ \sigma(r) = \frac{\rho\omega^2}{2} (L^2 - r^2) \qquad (C.26) $$

$$ u(r) = \frac{\rho r \omega^2}{2E} \left(L^2 - \frac{r^2}{3}\right) \qquad (C.27) $$

3. On s'intéresse à présent à la solution élastoplastique de ce problème. On définit un trajet de chargement monotone croissant pour la vitesse de rotation de la pâle. Comparer vos résultats à la solution numérique de CAST3M en traçant :

— le déplacement radial et la contrainte normale le long du rayon de la pâle *au dernier incrément de charge*,

— la courbe de charge de la structure. On tracera dans un graphe la vitesse de rotation de la pâle en fonction du déplacement de son extrémité. Que constatez-vous ?

4. En soustrayant les champs de contrainte calculés avec les solutions élastique et élastoplastique, montrer que la plasticité limite la croissance des contraintes dans la structure. On tracera ces résultats en fonction du rayon de la pâle.

5. On étudie à présent la décharge de la structure. On suppose que le turboréacteur s'arrête suffisamment doucement en passant par une succession d'état d'équilibre. Faire une reprise de calcul à partir du dernier pas de temps.

— La décharge est-elle élastique ?

— Tracer les champs résiduels (déplacement radial et la contrainte normale) le long du rayon de la pâle. Commenter.

## C.2 Système bidimensionnel axisymétrique

### C.2.1 Développement du code de calcul par éléments finis

**Formulation du problème éléments finis 2D axisymétrique et algorithme**

L'hypothèse 2D axisymétrique revient à définir un champ de déplacement de la forme :

$$ \mathbf{u}(r, z) = u_r(r, z)\mathbf{e}_r + u_z(r, z)\mathbf{e}_z \qquad (C.28) $$

où l'invariance du problème (géométrie, chargement, conditions aux limites et symétries matérielles) par rapport à la coordonnées $\theta$ a été prise en compte, il est en plus supposé
<figure><img src="image_83.png" /></figure>

qu'aucun chargement de torsion n'intervient dans l'étude ($u_\theta = 0$). Ceci conduit à des champs de déformation et de contrainte ayant quatre composantes non nulles :

$$ \sigma = \begin{cases} \sigma_{rr} \\ \sigma_{zz} \\ \sigma_{\theta\theta} \\ \sigma_{rz} \end{cases}; \quad \varepsilon = \begin{cases} \varepsilon_{rr} \\ \varepsilon_{zz} \\ \varepsilon_{\theta\theta} \\ 2\varepsilon_{rz} \end{cases} \qquad (C.29) $$

La résolution du problème d'évolution se fait toujours ici dans le cadre des petites perturbations. À chaque incrément de charge, le problème à résoudre prend toujours la forme (C.1); les forces externes et internes s'écrivent à présent :

$$ \mathbf{f}^{\text{ext}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{N}^T \mathbf{f}_d \, d\Omega + \int_{\partial\Omega_e \cap \partial\Omega_F} \mathbf{N}^T \mathbf{F}_d \, dS \right) \qquad (C.30) $$

$$ \mathbf{f}^{\text{int}} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{B}^T \sigma \, d\Omega \right) \qquad (C.31) $$

où la matrice **B** reliant les déformations aux déplacements nodaux est constituée par assemblage des matrices **B**<sub>p</sub> associées à chaque nœud **p** telle que :

$$ \mathbf{B}_p = \begin{bmatrix} N_{p,r} & 0 \\ 0 & N_{p,z} \\ \frac{N_p}{r} & 0 \\ N_{p,z} & N_{p,r} \end{bmatrix} \qquad (C.32) $$

À chaque itération *k* du calcul, le système linéaire suivant est résolu :

$$ \mathbf{K}^{(k)} \delta \mathbf{u}^{(k)} = \mathbf{R}^{(k)} \\ \text{avec} \quad \left| \mathbf{K}^{(k)} = - \frac{\partial \mathbf{R}}{\partial \mathbf{u}} \right|^{(k)} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \mathbf{B}^T \mathbf{H}^{(k)} \mathbf{B} \, d\Omega \right) \qquad (C.33) \\ \left| \mathbf{R}^{(k)} = \mathbf{f}^{\text{ext}} - (\mathbf{f}^{\text{int}})^{(k)} \right. $$

où $\mathbf{K}^{(k)}$ est la matrice tangente calculée à l'itération *k* du calcul, $\mathbf{H}^{(k)}$ l'opérateur tangent calculé lors de l'intégration de la loi de comportement élastoplastique. Toute intégrale de volume est effectuée avec les coordonnées cylindriques :

$$ \int_{\Omega} \dots d\Omega = 2\pi \int dz \int r dr \qquad (C.34) $$

L'algorithme C.1 reste inchangé, seuls les calcul et assemblage des matrices et des vecteurs forces doivent tenir compte de la dimension de l'espace. Des outils particuliers développés dans MatLab sont mis à votre disposition à cet effet.

**Outils MatLab pour la rédaction du code éléments finis**

Sont mises à votre disposition deux routines MatLab Matrix2D.m et Vector2D.m afin de faciliter le calcul des premiers et seconds membres en 2D. Elles permettent respectivement de calculer à partir d'un maillage bidimensionnel une matrice associée à une forme bilinéaire $a(u^*, u)$ et un vecteur associé à une forme linéaire $l(u^*)$. Une boucle sur les éléments, l'intégration des quantités élémentaires et leur assemblage sont réalisés dans ces routines. Les matrice et vecteur associés à ces formes bilinéaire et linéaire font intervenir les dérivées des fonctions de forme N par rapport aux
première et seconde coordonnées, ainsi que des champs définis aux nœuds (nodeField), constants par éléments (elemField), des champs définis comme une fonction arbitraire d'un champ défini aux nœuds (fctNodeField) et des champs définis aux points d'intégration des éléments (integPointField).

$$a(u^*, u) = \int_{\Omega} \frac{d^{(p)}u^*}{dx_i^{(p)}} \left( \prod_k \text{Fields} \right) \frac{d^{(q)}u}{dx_j^{(q)}} d\Omega \approx u^*Mu \quad (C.35)$$

$$l(u^*) = \int_{\Omega} \frac{d^{(p)}u^*}{dx_i^{(p)}} \left( \prod_k \text{Fields} \right) d\Omega \approx u^*f \quad (C.36)$$

où $(i, j) = (1, 2)$, $(x_1, x_2) \equiv (x, y)$ de sorte que la dérivée de chaque argument des formes bilinéaire et linéaire peut être prise par rapport aux coordonnées $x$ ou $y$. La routine `Matrix2D.m` prend en entrée les arguments suivants :

- `Coord`: table de coordonnées des nœuds du maillage, de taille $(N_{\text{nœuds}}, 2)$.

- `CONEC`: table de connectivité du maillage, de taille $(N_{\text{éléments}}, X)$, où $X$ est le nombre de nœuds de chaque élément.

- `op1x`: ordre de dérivation du premier argument $(u^*)$ de la forme bilinéaire par rapport à la première coordonnée $x$. 0 (rien) 1 (dérivation par rapport à $x$).

- `op1y`: comme `op1x` pour $y$.

- `op2x`: ordre de dérivation du second argument $(u)$ de la forme bilinéaire par rapport à la première coordonnée $x$. 0 (rien) 1 (dérivation par rapport à $x$).

- `op2y`: comme `op2x` pour $y$.

- `nodeField`: champs défini aux nœuds intervenant dans le calcul de l'intégrale. 1 champ par colonne, de taille $(N_{\text{nœuds}}, N_{\text{champs}})$.

- `elemField`: champs défini par éléments intervenant dans le calcul de l'intégrale. 1 champ par colonne, de taille $(N_{\text{éléments}}, N_{\text{champs}})$.

- `fctNodeField`: champs défini comme une fonction arbitraire d'un champ défini aux nœuds intervenant dans le calcul de l'intégrale. Si on a $N$ champs, on a une cellule de taille $N \times 2$ pour lequel chaque élément de la première ligne correspond aux valeurs nodales des champs, et la deuxième ligne est la fonction modifiant chaque champ (function handle)

- `integPointField`: champs défini aux points d'intégration d'un élément intervenant dans le calcul de l'intégrale de taille $(N_{\text{éléments}}, N_{\text{pointsIntégration}})$.

- `eltType`: type d'éléments finis utilisés : 'tri'(Triangle à 3 nœuds), 'quad'(quadrilatère à 4 nœuds).

Il est à noter que `op1x` et `op1y` sont complémentaires tout comme `op2x` et `op2y`. Si `op1x=1`, `op1y` doit être égal à zéro, et vice-versa. Les deux peuvent être égal à zéro, on définit alors une matrice de masse. La syntaxe est la suivante :

$$ 
\begin{aligned}
\mathbf{M} &= \text{Matrix2D}(\text{Coord}, \text{CONEC}, \text{op1x}, \text{op1y}, \text{op2x}, \text{op2y}, \text{nodeField}, \text{elemField}, \\
& \text{fctNodeField}, \text{integPointField}, \text{eltType}) \\
&= \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \left( \frac{d^{(\text{op1})}\mathbf{N}}{dx^{(\text{op1})}} \left( \prod_k \text{Fields} \right) \frac{d^{(\text{op2})}\mathbf{N}^T}{dx^{(\text{op2})}} \right) d\Omega \right)
\end{aligned}
\quad (C.37) 
$$
<figure><img src="image_84.png" /></figure>

La forme linéaire $l(u^*)$ n'ayant qu'un argument, la routine `Vector2D.m` ne nécessite de définir que les indices `opx` et `opy`. La syntaxe pour cette routine s'écrit :

$$ 
\begin{aligned}
\mathbf{f} &= \text{Vector2D}(\text{Coord}, \text{CONEC}, \text{opx}, \text{opy}, \text{nodeField}, \text{elemField}, \\
&\text{fctNodeField}, \text{integPointField}, \text{eltType}) \\
&= \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \frac{d^{(\text{op})}\mathbf{N}^T}{dx^{(\text{op})}} \left( \prod_k \text{Fields} \right) d\Omega \right)
\end{aligned}
\quad (C.38)
$$

Cependant, ces routines permettent de calculer les matrice et vecteur associés aux formes bilinéaire et linéaire ayant pour arguments des champs scalaires. Aussi pour la mécanique, il faut utiliser ces outils sur chacun des termes associés au développement des équations (C.29)-(C.33). Le calcul du vecteur des forces internes (C.31) provient d'une forme linéaire qui s'identifie ici à l'opposé du travail virtuel des forces intérieures :

$$ -W_{\text{int}}^* = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \sigma_{rr} \varepsilon_{rr}^* + \sigma_{zz} \varepsilon_{zz}^* + \sigma_{\theta\theta} \varepsilon_{\theta\theta}^* + 2\sigma_{rz} \varepsilon_{rz}^* \, d\Omega \right) \quad (C.39) $$

$$ = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \left( \sigma_{rr} u_{r,r}^* + \sigma_{zz} u_{z,z}^* + \sigma_{\theta\theta} \frac{u_r^*}{r} + \sigma_{rz} (u_{r,z}^* + u_{z,r}^*) \right) 2\pi r dr dz \right) \quad (C.40) $$

$$ = \sum_{e=1}^{N_e} \left( \sum_{p=1}^{N_n} \left\langle u_r^{(p)} \ u_z^{(p)} \right\rangle^* \int_{\Omega_e} \begin{bmatrix} N_{p,r} & 0 & \frac{N_p}{r} & N_{p,z} \\ 0 & N_{p,z} & 0 & N_{p,r} \end{bmatrix} \begin{bmatrix} \sigma_{rr} \\ \sigma_{zz} \\ \sigma_{\theta\theta} \\ \sigma_{rz} \end{bmatrix} 2\pi r dr dz \right) \quad (C.41) $$

$$ = \sum_{e=1}^{N_e} \left( \sum_{p=1}^{N_n} \left( \mathbf{u}^{(p)*} \right)^T \left( \mathbf{f}_p^{\text{int}} \right)^e \right) = \left( \mathbf{u}^{*} \right)^T \mathbf{f}_p^{\text{int}} \quad (C.42) $$

Pour vous faciliter la vie, la routine `InternalForces2Daxis.m` a été rédigée et utilise la routine `Matrix2D.m` en l'appliquant à l'équation (C.40). La syntaxe est la suivante :

$$ \mathbf{f}_{\text{int}} = \text{InternalForces2Daxis}(\text{Coord}, \text{CONECT}, \text{sigma}, \text{eltType}) \quad (C.43) $$

où `Coord`, `CONECT`, `sigma` et `eltType` sont les arguments de cette routine référant aux tables des coordonnées des nœuds et de connectivité du maillage, le champ de contrainte définis aux points d'intégration (`integPointField`) et le type d'élément fini utilisé respectivement.

D'autre part, le calcul de la matrice (élastique ou tangente (C.33)) nécessite aussi un développement particulier à partir des outils présentés ci-dessus. Étant donné un opérateur de rigidité (élastique ou tangent) **H**, cette matrice est telle que :

$$ (\mathbf{u}^*)^T \mathbf{Ku} = \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \langle \varepsilon_{rr} \ \varepsilon_{zz} \ \varepsilon_{\theta\theta} \ 2\varepsilon_{rz} \rangle^* \begin{bmatrix} H_{11} & H_{12} & H_{13} & H_{14} \\ H_{22} & H_{23} & H_{24} & H_{34} \\ \vdots & \vdots & \vdots & \vdots \\ H_{44} & H_{45} & H_{46} & H_{47} \end{bmatrix} \begin{Bmatrix} \varepsilon_{rr} \\ \varepsilon_{zz} \\ \varepsilon_{\theta\theta} \\ 2\varepsilon_{rz} \end{Bmatrix} 2\pi r dr dz \right) $$
$$
\begin{align*}
&= \sum_{e=1}^{N_e} \left( \int_{\Omega_e} \left[ \begin{aligned}[t]
& u_{r,r}^* (H_{11}u_{r,r} + H_{13}\frac{u_r}{r} + H_{14}u_{r,z}) + \frac{u_r^*}{r} (H_{31}u_{r,r} + H_{33}\frac{u_r}{r} + H_{34}u_{r,z}) \\
& + u_{r,z}^* (H_{41}u_{r,r} + H_{43}\frac{u_r}{r} + H_{44}u_{r,z}) \end{aligned} \right] 2\pi r dr dz \right. \\
&\qquad \left. + \int_{\Omega_e} \left[ u_{z,z}^* (H_{22}u_{z,z} + H_{24}u_{z,r}) + u_{z,r}^* (H_{42}u_{z,z} + H_{44}u_{z,r}) \right] 2\pi r dr dz \right. \\
&\qquad \left. + \int_{\Omega_e} \left[ u_{r,r}^* (H_{12}u_{z,z} + H_{14}u_{z,r}) + \frac{u_r^*}{r} (H_{32}u_{z,z} + H_{34}u_{z,r}) + u_{r,z}^* (H_{42}u_{z,z} + H_{44}u_{z,r}) \right] 2\pi r dr dz \right. \\
&\qquad \left. + \int_{\Omega_e} \left[ u_{z,z}^* (H_{21}u_{r,r} + H_{23}\frac{u_r}{r} + H_{24}u_{r,z}) + u_{z,r}^* (H_{41}u_{r,r} + H_{43}\frac{u_r}{r} + H_{44}u_{r,z}) \right] 2\pi r dr dz \right) \\
&= (\mathbf{u}_r^*)^T \mathbf{K}_{rr} \mathbf{u}_r + (\mathbf{u}_z^*)^T \mathbf{K}_{zz} \mathbf{u}_z + (\mathbf{u}_r^*)^T \mathbf{K}_{rz} \mathbf{u}_z + (\mathbf{u}_z^*)^T \mathbf{K}_{zr} \mathbf{u}_r \tag{C.44}
\end{align*}
$$

où les matrices $\mathbf{K}_{rr}$, $\mathbf{K}_{zz}$ et $\mathbf{K}_{rz} = \mathbf{K}_{zr}$ regroupent les termes de la matrice globale $\mathbf{K}$ associés respectivement aux équations $rr$, $zz$ et $rz$ de chaque nœud $p$ du maillage.

La routine `StiffnessMatrix2Daxis.m` a été rédigé et codée en utilisant la routine `Matrix2D.m` pour chacun des termes de l'équation (C.44). La syntaxe est la suivante :

$$ \mathbf{K} = \text{StiffnessMatrix2Daxis}(\text{Coord}, \text{CONECT}, \text{H}, \text{eltType}) \quad (\text{C.45}) $$

où les mêmes arguments que pour la routine `InternalForces2Daxis.m` sont requis, l'opérateur de rigidité $\mathbf{H}$ vient remplacer le champ de contrainte. Il doit être défini par une cellule contenant les 10 composantes indépendantes de la façon suivante :

$$ H^T = \{H_{11} \ H_{12} \ H_{22} \ H_{13} \ H_{23} \ H_{33} \ H_{14} \ H_{24} \ H_{34} \ H_{44}\} \quad (\text{C.46}) $$

Chaque composante est définie comme un `integPointField`, c'est-à-dire un tableau de taille ($N_{\text{éléments}}$, $N_{\text{pointsIntégration}}$).

Enfin, la routine `SurfaceTractions.m` permet de calculer le vecteur forces extérieures dans le cas 2D général, défini par le second terme de l'équation (C.30), associé à des forces surfaciques appliquées sur une partie du bord du domaine. La syntaxe de la routine est la suivante :

$$ \mathbf{f}_{\text{ext}} = \text{SurfaceTractions}(\text{Coord}, \text{CONECT}, \text{forces}) \quad (\text{C.47}) $$

où `Coord`, `CONECT`, `forces` sont respectivement la table des coordonnées des nœuds appartenant au bord sur lequel les forces surfaciques sont imposées, la table de connectivité de ce maillage 1D, et les forces surfaciques définies comme un `NodeField`, c'est-à-dire ici un tableau de taille ($N_{\text{éléments}}$, $N_{\text{pointsIntégration}}$). Attention cette routine doit être adaptée au cas axisymétrique, en adaptant la valeur du jacobien de la transformation qu'on pourra multiplier par la composante radiale de la force extérieure appliquée.

### C.2.2 Application : sphère sous pression interne

On considère un récipient sphérique creux (figure C.2) de rayon intérieur $a$, de rayon extérieur $b$, soumis à une pression interne $p$ croissante à partir de zéro. Le matériau considéré est élastoplastique et obéit au critère de von Misès. Il est isotrope élastiquement et plastiquement. Des écrouissages linéaires cinématique et isotrope seront successivement considérés dans l'analyse suivante. On exploite les symétries du problème (géométriques, de chargement et matérielles), aussi l'étude numérique sera conduite sur un quart de section de la sphère en option axisymétrique. La figure C.3 présente un premier maillage grossier du domaine d'étude réalisé avec le code CAST3M.
<figure><img src="image_85.png" /><figcaption>Figure C.2 – Sphère sous pression</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Figure C.3 – Maillage réglé grossier du domaine considéré</figcaption></figure>

**Partie 1 : sphère élastoplastique avec un écrouissage cinématique linéaire** Le fichier `sphereCinematique.dgibi` de la mise en données CAST3M de ce problème vous est fourni. Il sert ici de base de comparaison pour vos résultats. D'autre part, la solution analytique de cette sphère vous est aussi donnée pour un matériau élastique parfaitement plastique. Celle-ci permettra d'étudier l'influence de l'écrouissage sur la réponse élastoplastique de cette structure.

1. Faire le maillage du domaine de calcul. On prendra $a = 0.1\text{m}$ et $b = 0.2\text{m}$. Proposer plusieurs maillages et commenter le(s) type(s) d'élément(s) fini(s) choisi(s) et la zone de raffinement du maillage.

2. On considère une pression de 5000bars à l'intérieur du récipient sphérique, supposé en acier de construction classique S450 ($E = 210000\text{MPa}$, $\nu = 0.3$, $\sigma_0 = 450\text{MPa}$) ayant un module d'écrouissage (constant) de $h = 5000\text{MPa}$. Comparer les solutions numériques (CAST3M et la votre) et analytique purement élastiques. On rappelle ci-dessous la solution analytique purement élastique d'une sphère sous pression interne :

$$ 
\begin{align*}
\sigma_{rr} &= p \frac{a^3}{b^3 - a^3} \left(1 - \frac{b^3}{r^3}\right) \\
\sigma_{\theta\theta} &= \sigma_{\varphi\varphi} = p \frac{a^3}{b^3 - a^3} \left(1 + \frac{b^3}{2r^3}\right) \tag{C.48} \\
u_r &= \frac{p}{E} \frac{a^3}{b^3 - a^3} \left[(1 - 2\nu)r + (1 + \nu)\frac{b^3}{2r^2}\right]
\end{align*}
$$

3. On s'intéresse à présent à la réponse élastoplastique du système. On définit un trajet de chargement monotone croissant pour la pression appliquée à l'intérieur de la sphère.

(a) Comparer les solutions numériques (CAST3M et la votre) en traçant :

i. le déplacement radial et les composantes de contrainte radiale et orthoradiale le long du rayon de la pâle. Les composantes de contraintes et leur dérivée par rapport au rayon $r$ sont-elles continues ? Expliquez pourquoi.

ii. l'évolution du rayon plastique $c$, rayon correspondant à l'interface entre les zones élastique et élastoplastique, en fonction de la pression intérieure. Montrer que cette évolution est bien monotone croissante.

iii. la courbe de charge de la structure. On tracera dans un graphe l'évolution de la pression appliquée à l'intérieur de la sphère en fonction du déplacement radial d'un point de la paroi interne de la sphère. Que constatez-vous ?

(b) Comparer les solutions numériques (CAST3M et la votre) à la solution analytique donnée pour un matériau élastique parfaitement plastique. On rappelle l'expression
analytique du champ de contrainte dans la zone plastique :

$$
\begin{align}
\sigma_{rr} &= 2\sigma_0 \ln\left(\frac{r}{c}\right) - \frac{2}{3}\sigma_0 \left(1 - \frac{c^3}{b^3}\right) \tag{C.49} \\
\sigma_{\theta\theta} &= \sigma_{\varphi\varphi} = \sigma_{rr} + \sigma_0 \nonumber
\end{align}
$$

où c est le rayon de la zone élastoplastique, solution de l'équation suivante :

$$
p = 2\sigma_0 \ln\left(\frac{c}{a}\right) + \frac{2}{3}\sigma_0 \left(1 - \frac{c^3}{b^3}\right) \quad (C.50)
$$

et celle du champ de contrainte dans la zone élastique :

$$
\begin{align}
\sigma_{rr} &= \frac{2}{3}\sigma_0 \frac{c^3}{b^3} \left(1 - \frac{b^3}{r^3}\right) \tag{C.51} \\
\sigma_{\theta\theta} &= \sigma_{\varphi\varphi} = \frac{2}{3}\sigma_0 \frac{c^3}{b^3} \left(1 + \frac{b^3}{2r^3}\right) \nonumber
\end{align}
$$

i. Dans un graphe, superposez les solutions des composantes de contrainte radiale et orthoradiale données numériquement avec écrouissage et analytiquement sans écrouissage, le long du rayon de la pâle. Quel pourcentage d'erreur y a-t-il sur les contraintes maximales entre les solutions obtenues avec et sans écrouissage. Que peut-on en conclure ?

ii. Faire la même chose avec le déplacement radial.

**Attention** L'option axisymétrique avec laquelle vous travaillez vous donne les contraintes dans le système de coordonnées cylindriques ; les comparaisons doivent être effectuées sur les composantes de contraintes exprimées dans le système de coordonnées sphériques (voir annexe C.3).

4. En soustrayant les champs de contrainte calculés avec les solutions élastique et élastoplastique (avec CAST3M et votre code), montrer que la plasticité limite la croissance des contraintes dans la structure. On tracera ces résultats en fonction du rayon de la sphère.

5. On étudie à présent la décharge de la structure. Faire une reprise de calcul à partir du dernier pas de temps. La décharge est-elle élastique ? Superposez sur un graphe les composantes de contraintes résiduelles obtenues avec CAST3M et vos résultats.

6. Les contraintes résiduelles données par la solution analytique (obtenues pour un *matériau élastique parfaitement plastique*) sont les suivantes :

$$
\begin{align}
\sigma_{rr} &= 2\sigma_0 \left[ \ln \frac{r}{b} - \frac{a^3}{b^3 - a^3} \left(1 - \frac{b^3}{r^3}\right) \ln \frac{b}{a} \right] \tag{C.52} \\
\sigma_{\theta\theta} &= \sigma_{\varphi\varphi} = \sigma_0 \left[ 2 \ln \frac{r}{b} + 1 - 2 \frac{a^3}{b^3 - a^3} \left(1 + \frac{b^3}{2r^3}\right) \ln \frac{b}{a} \right] \nonumber
\end{align}
$$

La prise en compte de l'écrouissage change-t-il beaucoup l'état résiduel prédit ? Superpo-
sez les champs sur un même graphe et comparer les valeurs maximales.

**Partie 2 : sphère élastoplastique à écrouissage isotrope linéaire** Le fichier de la mise en données CAST3M de ce problème est fourni : sphereIsotrope.dgibi. Il prend en compte un écrouissage isotrope linéaire.

1. Refaire la même étude avec un matériau à écrouissage isotrope.

2. Comparer les états résiduels (contraintes et distorsions) obtenus respectivement avec les deux types d'écrouissage. Commenter.
<figure><img src="image_87.png" /></figure>

C.3 Composantes de contrainte en coordonnées sphériques

Le calcul mené en option 2D axisymétrique vous donne les champs exprimés dans la base cylindrique (e_r, e_θ, e_z). Il vous faut les exprimer dans le système de coordonnées sphériques pour effectuer la comparaison avec la solution analytique d'une part, et la solution numérique obtenue avec le code CAST3M d'autre part.

Les composantes de contrainte radiale et orthoradiale dans la base sphérique s'expriment en fonction des composantes de contrainte dans la base cylindrique de la façon suivante : soit e'_r le vecteur unitaire lié à la base sphérique, tel que e'_r · e_r = cos φ, alors la composante radiale dans la base sphérique s'obtient comme suit :

$$
\begin{align}
\sigma'_{rr} &= \mathbf{e}'_r \cdot \boldsymbol{\sigma} \cdot \mathbf{e}'_r \nonumber \\
&= \langle \cos \varphi \ 0 \ \sin \varphi \rangle \begin{bmatrix} \sigma_{rr} & 0 & \sigma_{rz} \\ 0 & \sigma_{\theta\theta} & 0 \\ \sigma_{rz} & 0 & \sigma_{zz} \end{bmatrix}_{(e_r, e_\theta, e_z)} \begin{Bmatrix} \cos \varphi \\ 0 \\ \sin \varphi \end{Bmatrix} \tag{C.53} \\
&= \sigma_{rr} \cos^2 \varphi + \sigma_{zz} \sin^2 \varphi - 2\sigma_{rz} \cos \varphi \sin \varphi \nonumber
\end{align}
$$

La composante orthoradiale est quant elle obtenue de manière analogue :

$$
\begin{align}
\sigma_{\varphi\varphi} &= \mathbf{e}_{\varphi} \cdot \boldsymbol{\sigma} \cdot \mathbf{e}_{\varphi} \nonumber \\
&= \langle -\sin\varphi \ 0 \ \cos\varphi \rangle \begin{bmatrix} \sigma_{rr} & 0 & \sigma_{rz} \\ 0 & \sigma_{\theta\theta} & 0 \\ \sigma_{rz} & 0 & \sigma_{zz} \end{bmatrix}_{(e_r, e_\theta, e_z)} \begin{Bmatrix} -\sin\varphi \\ 0 \\ \cos\varphi \end{Bmatrix} \tag{C.54} \\
&= \sigma_{rr} \sin^2\varphi + \sigma_{zz} \cos^2\varphi - 2\sigma_{rz} \cos\varphi \sin\varphi \nonumber
\end{align}
$$
Bibliographie

[1] Pierre SUQUET. *Rupture et plasticité*. École Polytechnique, 2004 (cf. pages 5, 6, 40).

[2] Jean SALENÇON. *De l'élastoplasticité au calcul à la rupture*. Éditions de l'École poly-technique, 2002. ISBN : 9782730209151 (cf. pages 7, 19, 40).

[3] Jacques BESSON, Georges CAILLETAUD, Jean-Louis CHABOCHE et Samuel FOREST. *Mécanique non linéaire des matériaux*. Hermès, 2001. ISBN : 9782746202689 (cf. page 7).

[4] Jean LEMAITRE, Jean-Louis CHABOCHE, Ahmed BENALLAL et Rodrigue DESMORAT. *Mécanique des matériaux solides*. 3e édition. Dunod, 2009. ISBN : 9782100516230 (cf. page 7).

[5] Huy Duong BUI. « Évolution de la frontière du domaine élastique des métaux avec l'écrouissage plastique et le comportement élastoplastique d'un agrégat de cristaux cubiques ». *Mémoire de l'Artillerie Française* 1 (1971), pages 141-165 (cf. page 11).

[6] Rodney HILL. *The Mathematical Theory of Plasticity*. Clarendon Press, 1967. ISBN : 9780198503675 (cf. page 15).

[7] Jean-Jacques MARIGO. *Plasticité et rupture*. [cel-01374813]. École Polytechnique, 2011 (cf. pages 18, 40).

[8] Juan Carlos SIMO et Thomas HUGHES. *Computational Inelasticity*. Springer, 1998. ISBN : 9780387975207 (cf. page 67).

[9] Paul GERMAIN. *Mécanique. Cours X*. Ellipses, 1986. ISBN : 9782729886530.

[10] Patrick De BUHAN. *Plasticité et calcul à la rupture*. Presses Ponts et Chaussées, 2007. ISBN : 9782859784324.

[11] Olgierd ZIENKIEWICZ, Robert TAYLOR et David FOX. *The Finite Element Method for Solid and Structural Mechanics*. 7e édition. Butterworth-Heinemann, 2014. ISBN : 9780080951362.

[12] Ted BELYTSCHKO, Wing Kam LIU, Brian MORAN et Khalil ELKHODARY. *Nonlinear Finite Elements for Continua and Structures*. 2e édition. Wiley, 2000. ISBN : 9781118632703.

[13] Thomas HUGHES. *The Finite Element Method. Linear Static and Dynamic Finite Element Analysis*. Dover, 1987. ISBN : 9780486411811.

[14] Paul GERMAIN, Pierre SUQUET et Quoc Son NGUYEN. « Continuum thermodynamics ». *Journal of Applied Mechanics* 50 (1983), pages 1010-1020. DOI : 10.1115/1.3167184.

[15] Bernard JAOUL. *Étude de la plasticité et application aux métaux*. Mines ParisTech, 2008. ISBN : 9782911762918.

[16] Hervé OUDIN. *Introduction à la plasticité*. [cel-00483603]. École Centrale de Nantes, 2009.

[17] Pierre SUQUET. « Plasticité et homogénéisation ». Thèse de doctorat. Université Pierre et Marie Curie, 1982.
<figure><img src="image_88.png" /></figure>

[18] François SIDOROFF. *Comportement des matériaux*. École Centrale de Lyon, 1984.

[19] Jean-Baptiste LEBLOND. *Plasticité et viscoplasticité en petites et grandes transformations*. Université Paris VI, 2007.

[20] René BILLARDON. *Comportement des matériaux*. Cours de DEA. Université Paris VI, 2007.

[21] Jean MANDEL. *Plasticité classique et viscoplasticité*. Tome 97. CISM International Centre for Mechanical Sciences. Springer-Verlag, 1971. ISBN : 9783211811979.