# Comportement élastoviscoplastique d'une fonte GS et son intégration numérique

Fabien Szmytka, Habibou Maitournam, Luc Rémy, Alain Köster, Myriam Bourgeois, Thi Mac-Lan Nguyen-Tajan

► To cite this version:

Fabien Szmytka, Habibou Maitournam, Luc Rémy, Alain Köster, Myriam Bourgeois, et al. Comportement élastoviscoplastique d'une fonte GS et son intégration numérique. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. pp.283-288. hal-01487614

HAL Id: hal-01487614

https://hal.science/hal-01487614v1

Submitted on 12 Mar 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Comportement élastoviscoplastique d'une
fonte GS et son intégration numérique

Fabien Szmytka*,**,*** — Habibou Maitournam** — Luc Rémy***
— Alain Köster*** — Myriam Bourgeois* — T. Mac Lan Nguyen-
Tajan*

* PSA Peugeot Citroën - Direction de la Recherche et de l'Innovation Automobile
Route de Gisy - 78 943 Vélizy-Villacoublay Cedex

** Laboratoire de Mécanique des Solides (CNRS UMR 7649)
École Polytechnique - 91 128 Palaiseau Cedex

*** Centre des Matériaux P.M. Fourt (CNRS UMR 7633)
École Nationale Supérieure des Mines de Paris, BP87 - 91003 Evry Cedex
fabien.szmytka@mpsa.com, habibou@lms.polytechnique.fr, luc.remy@mat.ensmp.fr

RÉSUMÉ. Le travail présenté ici porte sur la détermination et l’intégration numérique d’une loi de comportement mécanique élasto-viscoplastique. Cette dernière est développée dans le but de dimensionner à la fatigue thermomécanique des structures de l’industrie automobile en décrivant au mieux la viscosité sur une large gamme de sollicitation. Ce nouveau modèle a été implémenté dans le code Abaqus en se basant sur l’algorithme du retour radial. Des simulations thermomécaniques ont été ensuite réalisées sur des collecteurs d’échappement et ont permis de statuer avec succès sur la durée de vie de la structure.

ABSTRACT. A new elasto-viscoplastic behaviour law and its numerical integration scheme have been developed. Our aim is to design hot automotive components with a good description of all the viscous effects for a wide load range. This new model has been integrated in the Abaqus code thanks to a return-mapping algorithm. Thermo-mechanical simulations performed on exhaust manifolds had led to good estimation of the component lifetime.

MOTS-CLÉS : loi de comportement élasto-viscoplastique anisotherme, intégration numérique, algorithme du retour radial, fonte à graphite sphéroïdal, fatigue oligocyclique, collecteur d'échappement

KEYWORDS: Elasto-viscoplastic anisotherm behaviour law, numerical integration, return mapping algorithm, spheroidal cast-iron, low-cycle fatigue, exhaust manifold
# 1. Introduction

La réduction de la taille des cylindrées des moteurs diesel pour diminuer la consommation de carburant ainsi que les diverses stratégies de dépollution soumettent des pièces comme la culasse ou le collecteur d'échappement à des cycles thermiques de plus en plus sévères pouvant entraîner des dérives dimensionnelles par accumulation de déformation anélastique. Le dimensionnement de ces structures métalliques, travaillant à hautes températures, fait de plus en plus souvent appel à des lois de comportement élasto-viscoplastique. Ces lois macroscopiques, parmi lesquelles on peut citer celles issues des travaux de Lemaitre et Chaboche (Lemaitre et al., 1996), doivent décrire avec précision l'évolution des déformations et des contraintes qui sont à la base des critères de dimensionnement. Néanmoins, ces modèles ne parviennent souvent pas à décrire efficacement le comportement visqueux sur une large gamme de température (de l'ambiante à 700°C) et de sollicitation mécanique.

Ainsi, la modélisation du comportement de la fonte GS des collecteurs d'échappement a été jusqu'alors abordée en séparant les déformations visqueuses des déformations plastiques et en s'attachant à bien décrire les zones de chargement maximale en température et contrainte où apparaissent les fissures de fatigue (Charkaluk, 1999). Cependant, une telle loi ne permet pas de connaître avec précision la répartition des déformations inélastiques en tout point du collecteur et les dérives dimensionnelles qu'elles occasionnent sur la structure globale. La mise au point de modèles de comportement prenant mieux en compte l'ensemble des phénomènes visqueux agissant sur les matériaux et permettant, après implémentation numérique, d'aboutir à des calculs de structure prédictifs et d'une durée compatible avec les exigences de bureaux d'étude industriels est donc une étape essentielle du dimensionnement de ces pièces.

Nous nous intéressons ici au développement d'une telle loi de comportement depuis la phase expérimentale qui nous permet d'envisager la modélisation à adopter jusqu'à l'écriture des équations du modèle et leur implémentation dans un code de calcul éléments finis, étape sur laquelle nous nous attarderons plus longuement et que nous illustrerons à travers l'étude d'un collecteur d'échappement.

# 2. Stratégie Expérimentale

La constitution d'une base expérimentale robuste nous permet de caractériser le matériau, de mettre au point une loi de comportement et d'en identifier les paramètres. Afin de nous permettre d'appréhender au mieux le comportement visqueux de la fonte GS considérée<fnref n="1" />, nous avons principalement réalisé des essais d'écrouissage cyclique. Ces derniers, isothermes et pilotés en déformation à fréquence constante, sont effectués à différentes températures et afin de couvrir une gamme de sollicitation représentative de la structure. Pour un niveau de déformation $ε_1$ donné, ils permettent ainsi (voir figure 1) de couvrir une très large gamme de vitesse de déformation

<footnote n="1">1. contenant principalement du silicium et du molybdène et dénommée par la suite SiMo</footnote>
inélastique (de $10^{-3}s^{-1}$ à $10^{-9}s^{-1}$).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Essai d'écrouissage cyclique</figcaption></figure>

Pour une température donnée, chacune des phases a alors une durée fixe quel que soit $\varepsilon_1$. On sollicite l'éprouvette d'abord pour des valeurs de $\varepsilon_1$ faibles de l'ordre de 0, 1% et une fois la boucle contraintes déformations stabilisée, on augmente cette valeur. Un essai complet, se composant d'environ 8 à 10 paliers, est ainsi très riche d'information et, complété par des essais de tractions simples, nous permet de bien appréhender le comportement viscoplastique de la fonte SiMo. Nous pouvons alors envisager sa modélisation, qui pourra être ensuite validée notamment par des essais anisothermes uniaxiaux (Rémy, 2003).

### 3. Vers un nouveau modèle de comportement

Nous nous plaçons ici dans le cadre de la viscoplasticité unifiée. Le matériau obéisant au critère de von Mises, nous avons recours à un écrouissage cinématique non linéaire selon le modèle d'Armstrong-Frederik (Armstrong et al., 1966) :

$$ \underline{X} = \frac{2C}{3}\underline{\alpha} \qquad [1] $$

$$ \underline{\alpha} = \dot{\underline{\varepsilon}}_{vp} - \frac{3\dot{\gamma}D}{2C}\underline{X} \qquad [2] $$

avec $C$ et $D$, 2 paramètres à identifier, $\underline{X}$ la variable d'écrouissage et $\underline{\alpha}$, sa variable associée. La base expérimentale met de plus en évidence un adoucissement cyclique de la fonte SiMo. Nous avons donc adjoint un écrouissage isotrope selon le modèle de Chaboche (Lemaitre et al., 1996) :

$$ \dot{R} = b(Q - R)\dot{\gamma}, \text{ avec } R(t = 0) = R_0 \qquad [3] $$

avec $b$, $Q$ et $R_0$, 3 paramètres supplémentaires et $\gamma$, le multiplicateur viscoplastique. Les phases de relaxation des essais d'écrouissage cyclique présentent une évolution quasi linéaire de la contrainte en fonction de la vitesse de déformation inélastique dans un diagramme semi-logarithmique. En nous basant sur les mécanismes de déformation à l'échelle microscopique (Poirier, 1976), nous avons mis au point un nouveau potentiel viscoplastique :

$$ \dot{\gamma} = P(J_2(\underline{\underline{\sigma}} - \underline{\underline{X}}) - R) \left( \sinh \left( \frac{J_2(\underline{\underline{\sigma}} - \underline{\underline{X}}) - R}{K} \right) \right) \qquad [4] $$
avec *K*, un nouveau paramètre. Ce potentiel s'écrit donc comme la combinaison d'un terme en sinus hyperbolique valant pour les plus hautes températures et qui a montré son efficacité pour décrire le fluage à température élevée ((Chateau *et al.*, 2001) d'après (Dyson *et al.*, 1998)) et d'un terme multiplicatif polynomial en $\underline{\underline{\sigma}}$, *P*, lié à la viscoplasticité à basse température. La loi de comportement ainsi définie rentre dans le cadre thermodynamique des matériaux standards généralisés (Halphen *et al.*, 1975).

## 4. Intégration Numérique : algorithme de retour radial

Une fois les paramètres de la loi identifiés (nous obtenons une bonne corrélation expérimentale à toutes les températures. Les résultats à 500°C sont représentés sur la figure 2), il convient d'écrire son schéma d'intégration numérique afin de réaliser des calculs avec le code éléments finis *Abaqus*.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Essai d'écrouissage cyclique à 500°C - les traits pleins représentent la simulation</figcaption></figure>

Nous utilisons alors un schéma implicite utilisant la technique du retour radial (Nguyen, 1977, Simo *et al.*, 1985). Après discrétisation des variables du problème, deux étapes présentent certaines difficultés. Il nous faut tout d'abord définir l'incrément du multiplicateur viscoplastique noté $\Delta\gamma$. La règle de normalité s'exprime à l'instant $t_{n+1}$ selon<fnref n="2" /> :

$$ \Delta\underline{\underline{\epsilon}}^{vp} = \Delta\gamma \frac{\underline{\underline{A}}_{n+1}}{\| \underline{\underline{A}}_{n+1} \|} = \Delta\gamma \underline{n} \quad [5] $$

<footnote n="2">2. Toutes les valeurs relatives à l'instant $t_n$, connu, sont indicées par *n* et respectivement toutes celles relatives à l'instant $t_{n+1}$, que nous souhaitons calculer, par $n + 1$. Un incrément entre ces deux instants est référencé par un $\Delta$.</footnote>
avec $\| \underline{A}_{n+1} \| = \sqrt{\underline{A}_{n+1} : \underline{A}_{n+1}}$. Ainsi, sachant que $\underline{A}_{n+1} = \text{dev}(\underline{\sigma}_{n+1}) - \underline{X}_{n+1}$ et en écrivant les incréments de contrainte et des variables d'écrouissage, on obtient :

$$ \| \underline{A}_{n+1} \| = \| \underline{B}^*_{n} \| - (2\mu_{n+1} + 2\Delta\mu)\Delta\gamma - \frac{\frac{2}{3}(\Delta C + C_{n+1})\Delta\gamma}{\left(1 + \sqrt{\frac{2}{3}}D_{n+1}\Delta\gamma\right)} \quad [6] $$

avec $\underline{B}^*_{n}$, une fonction entièrement déterminée par la connaissance de l'ensemble des variables au temps $t_n$ ; $\lambda$ et $\mu$, les coefficients de Lamé. En combinant les équations 4 et 6, nous obtenons une formulation implicite en $\Delta\gamma$ qui est résolue grâce à un schéma itératif de Newton-Raphson d'ordre 1 et pouvons ainsi notamment expliciter $\underline{\Delta\underline{\epsilon}}^{vp}$.

Le secont point délicat réside dans le calcul de l'opérateur tangent cohérent qui détermine les variations de $\underline{\Delta\underline{\sigma}}$ en fonction de celles de $\underline{\Delta\underline{\epsilon}}$. Cette étape se résume à expliciter la dérivée de $\underline{\Delta\gamma}$ par rapport à $\underline{\Delta\underline{\epsilon}}$. Pour y parvenir, on considère $\underline{\Delta\gamma}$ et $\| \underline{A}_{n+1} \|$ comme deux fonctions à plusieurs variables à savoir $\underline{\Delta\gamma} = \eta_1(\| \underline{A}_{n+1} \|, R_{n+1})$ et $\| \underline{A}_{n+1} \| = \eta_2(\underline{\Delta\underline{\epsilon}}, \underline{\Delta\gamma})$. De cette façon, en écrivant leur différentielle, nous pouvons alors obtenir

$$ \frac{d\underline{\Delta\gamma}}{d\underline{\Delta\underline{\epsilon}}} = \frac{\left| \frac{\partial\eta_1}{\partial\|\underline{A}_{n+1}\|} \right|_{R_{n+1}} \cdot \left| \frac{\partial\eta_2(\underline{\Delta\underline{\epsilon}}, \underline{\Delta\gamma})}{\partial\underline{\Delta\underline{\epsilon}}} \right|_{\underline{\Delta\underline{\gamma}}} \quad [7] \\ 1 - \left( 1 - \left| \frac{\partial\eta_1}{\partial\|\underline{A}_{n+1}\|} \right|_{R_{n+1}} \cdot \left| \frac{\partial\eta_2(\underline{\Delta\underline{\epsilon}}, \underline{\Delta\gamma})}{\partial\underline{\Delta\gamma}} \right|_{\underline{\Delta\underline{\epsilon}}} - \left| \frac{\partial\eta_1}{\partial R_{n+1}} \right|_{\|\underline{A}_{n+1}\|} \cdot \frac{\partial R_{n+1}}{\partial\underline{\Delta\gamma}} \right) $$

qui s'explicite facilement. Ainsi, nous parvenons à écrire simplement l'opérateur tangent cohérent et à compléter le développement d'une sous-routine matériau Umat dans *Abaqus*. Ceci nous permet d'intégrer la nouvelle loi de comportement dans des calculs élément finis.

## 5. Applications : calcul d'un collecteur d'échappement

Un calcul thermomécanique (Lederer *et al.*, 2000) est réalisé sous *Abaqus* pour 10 cycles de sollicitations suivi d'un retour à température ambiante. Les simulations effectuées (voir figure 3), dans des temps de calcul comparables avec les lois classiques, présentent une bien meilleure corrélation expérimentale en terme de déformée globale.

## 6. Conclusions et perspectives

Nous avons ici succinctement présenté les différentes étapes du développement d'une nouvelle loi de comportement élasto-viscoplastique depuis la création d'une base expérimentale jusqu'à l'écriture du schéma d'intégration numérique et le calcul élément fini. Ce modèle macroscopique, basé sur un potentiel visqueux innovant, nous
<figure><img src="image_4.png" /><figcaption>**Figure 3.** *Calcul EF thermomécanique d'un collecteur d'échappement - Déformée après retour à l'ambiante*</figcaption></figure>

permet aujourd'hui de fournir des calculs de structures donnant des résultats en terme de déformée globale plus conformes à la réalité que les lois classiquement utilisées.

**7. Bibliographie**

Armstrong P., Frederick C., A mathematical representation of the multiaxial Bauschinger effet., Technical report, Berkeley Nuclear Laboratories, 1966.

Charkaluk E., Dimensionnement des structures à la fatigue Thermomécanique, PhD thesis, École polytechnique, 1999.

Chateau E., Rémy L., « Oxidation-assisted Creep Damage in Udimet 720LI », in, M. Associazione Italiana di Metallurgica (ed.), *EUROMAT 2001, proceedings of the 7th European Conference on Advanced Materials and Processes, RIMini, Italy*, 2001.

Dyson B., McLean M., *Microstructural stability of creep resistant alloys for high temperature plant applications*, A. Strang, J. Cawley, G.W. Greenwood Eds, Institute of Materials, Londres, chapter Microstructural evolution and its effects on the creep performance of high temperature alloys, p. 371-392, 1998.

Halphen B., Nguyen Q., « Sur les matériaux standards généralisés », *Journal de Mécanique*, vol. 14, n° 1, p. 39-63, 1975.

Lederer G., Charkaluk E., Verger L., Constantinescu A., « Numerical lifetime assesment of engine parts submitted to thermomecanical fatigue », *SAE Technical paper series*, 2000.

Lemaitre J., Chaboche J.-L., *Mécanique des matériaux solides*, 2ème edn, Dunod, 1996.

Nguyen Q., « On the elastic plastic initial-value problem and its numerical integration », *International Journal for Numerical Method in Engineering*, 1977.

Poirier J.-P., *Plasticité à haute température des solides cristallins*, Édition Eyrolles, 1976.

Rémy L., *Comprehensive Structural Integrity*, vol. 5, Elsevier, chapter Thermal-Mechanical Fatigue, p. 113-200, 2003.

Simo J., Taylor R., « Consistent tangent operators for rate-independent elastoplasticity », *Computer methods in applied mechanics and engineering*, vol. 48, p. 101-118, 1985.