# Simulation numérique du refendage de tôles minces

Nabil Belamri, Abel Cherouat, Khemais Saanouni, Philippe Autesserre

► To cite this version:

Nabil Belamri, Abel Cherouat, Khemais Saanouni, Philippe Autesserre. Simulation numérique du refendage de tôles minces. 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. hal-02572290

HAL Id: hal-02572290

https://utt.hal.science/hal-02572290

Submitted on 26 Apr 2023

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>

Distributed under a Creative Commons Attribution 4.0 International License
Simulation numérique 3D
du procédé de refundage de tôles minces

N. Belamri*, A. Cherouat*, K. Saanouni* et P. Autesserre**

* LASMIS-CNRS-FRE 2719, Université de technologie de Troyes
BP 2060 - 10010 Troyes - France,

** Centre Opérations Industrielles, Arcelor Research SA
Voie Romaine BP30320 57283 Maizières-lès-Metz

Résumé :

Nous présentons ici un modèle élastoplastique anisotrope intégrant l'écrouissage isotrope et cinématique non linéaire, fortement couplé à l'endommagement ductile isotrope. Ce modèle est implanté dans le code éléments finis Abaqus®. Les aspects théoriques et numériques du modèle sont présentés, un solveur dynamique explicite de résolution est notamment utilisé pour résoudre le problème d'équilibre. Une attention spéciale a été apportée à l'intégration des équations constitutives couplées, qui s'effectue à l'aide d'un schéma itératif implicite. Les capacités du modèle sont illustrées sur la simulation 3D du procédé de refundage d'acier. Le refundage est un procédé de découpe longitudinale en continu de tôles minces, le plus souvent installé dans les centres de service acier.

Abstract :

We present here an anisotropic elastoplastic model accounting for non-linear isotropic and kinematic hardening, strongly coupled
with isotropic ductile damage. This model is integrated in finite elements code Abaqus®. Both theoretical and numerical aspects
are introduced, an explicit dynamic resolution strategy is especially used to solve the associated initial and boundary value
problem. Moreover, a special care is given to the time integration of the fully coupled constitutive equations thanks to an implicit
iterative scheme. We illustrate the capabilities of the model on the 3D simulation of slitting. Slitting is a continuous metal
cutting process, usually installed in the steel service centers.

Mots clefs :

Refendage, anisotropie de plasticité, Critère de Hill, Endommagement ductile, Intégration implicite, Résolution dynamique.

1 Introduction

Un modèle de comportement élastoplastique anisotrope (critère de Hill) couplé à un endommagement ductile isotrope a été développé au LASMIS [1, 2]. Ce modèle décrit les principaux phénomènes physiques (écrouissage, endommagement, rupture) caractérisant le comportement mécanique des pièces soumises à de grandes déformations. L'importance du couplage fort sur la représentativité du modèle (comparaisons expérimentales) a été montré.

Après une présentation du modèle, nous décrivons une des applications possibles dans le domaine de la découpe des aciers. Nous simulons ainsi en 3D le procédé industriel de refundage de tôles [1].

Le refundage est un procédé de découpe longitudinale de tôle à l'aide de couteaux circulaires, permettant d'obtenir des bobines de largeur finale inférieure à la largeur de la bobine initiale. La machine utilisée est appelée une refendeuse, sa partie travaillante est constituée de deux arbres cylindriques (tête de refundage) sur lesquels sont positionnés les couteaux de découpe (Figure 1). Le maintien de ces couteaux est assuré par des entretoises, dont certaines sont revêtues de néoprène. Ces dernières, alors appelées entretoises dévêtisseuses, ont aussi pour fonction d'assurer le maintien de la tôle, tout en permettant une meilleure décohésion au moment de la coupe [1].
2 Formulations mécaniques et numériques du modèle

2.1 Equations constitutives

Le comportement élastoplastique couplé à l'endommagement est décrit dans le cadre de la
thermodynamique des processus irréversibles avec variables d'état [2]. Nous limitons la formulation au cas
isotherme, avec l'hypothèse d'une seule surface de charge pour la plasticité et l'endommagement isotrope.

Les phénomènes physiques pris en compte peuvent être représentés par des variables d'état (observables ou internes) définies comme suit :

* $(\underline{D}, \underline{\sigma})$ représentent le tenseur de taux de déformation d'Euler $\underline{D}$ associé au tenseur de contraintes de Cauchy $\underline{\sigma}$. Ce dernier est déduit de l'élasticité linéaire isotrope définie par l'équation (1). En petites déformations élastiques, le taux de déformation totale $\underline{D}$ (partie symétrique du gradient de vitesse de déformation) est décomposé en $\underline{D} = \dot{\underline{\epsilon}}^{ej} + \underline{D}^p$ sachant que $\dot{\underline{\epsilon}}^{ej}$ est le taux de déformation élastique de Jaumann et $\underline{D}^p$ est la partie plastique donnée par l'équation (6),

• $(\dot{\underline{\alpha}}, \underline{\underline{X}})$ variables tensorielles représentant l'écrouissage cinématique définies par les équations (7) et (2),

• $(\dot{r}, R)$ variables scalaires représentant l'écrouissage isotrope définies par les équations (8) et (2),

* $(\dot{\underline{D}}, \underline{Y})$ variables scalaires représentant l'endommagement ductile et la force thermodynamique associée définies par les équations (9) et (4).

a. variables de contrainte

$$
\underline{\sigma} = (1 - D)[\lambda \operatorname{tr}(\underline{\epsilon}^e)I + 2\mu \underline{\epsilon}^e] = (1 - D)\underline{\Lambda}: \underline{\epsilon}^e \quad (1)
$$

$$
\underline{X} = \frac{2}{3} C(1 - D)\underline{\alpha} \qquad (2)
$$

$$
R = Q(1 - D)r \tag{3}
$$

$$
Y = \frac{1}{2} \underline{\epsilon}^e : \underline{\Lambda} : \underline{\epsilon}^e + \frac{1}{3} C \underline{\alpha} : \underline{\alpha} + \frac{1}{2} Q r^2 \quad (4)
$$

b. Fonction de charge anisotrope

$$
f(\underline{\sigma}, \underline{\underline{X}}, R) = \frac{1}{\sqrt{1-D}} \|\underline{\sigma} - \underline{\underline{X}}\|_{H} - \frac{R}{\sqrt{1-D}} - \sigma_y \le 0 \quad \text{Sachant que } \|\underline{\sigma} - \underline{\underline{X}}\|_{H} = \sqrt{(\underline{\sigma} - \underline{\underline{X}}) : H : (\underline{\sigma} - \underline{\underline{X}})} \quad (5)
$$

c. *Variables de déformation*

$$
\underline{D}^p = \dot{\underline{\delta}}\underline{n}, \text{ avec } \underline{n} = \frac{\tilde{n}}{\sqrt{1-D}} = \frac{\underline{H} : (\underline{\sigma} - \underline{\underline{X}})}{\sqrt{1-D} \|\underline{\sigma} - \underline{\underline{X}}\|} \quad (6)
$$

$$
\dot{\underline{\alpha}} = \underline{D}^p - a \dot{\underline{\delta}} \underline{\alpha} \qquad (7)
$$

$$
\dot{r} = \dot{\delta}(\frac{1}{\sqrt{(1-D)}} - br) \qquad (8)
$$

$$
\dot{\underline{D}} = \dot{\delta}(\frac{Y - Y_0}{A})^{\gamma} \frac{1}{(1 - D)^{\beta}} \qquad (9)
$$

<caption>Tableau 1 : Présentation des équations constitutives couplées</caption>

Le Tableau 1 récapitule les équations constitutives couplées dérivant d'un potentiel dissipatif d'état [1, 2, 3]. Dans ces équations, on retrouve les constantes classiques de Lamé définies dans l'élasticité isotrope par λ et μ, le module d'écrouissage cinématique C, le module d'écrouissage isotrope Q, les coefficients non-linéaires
définissant la vitesse à laquelle le matériau s'écrouit a et b, le module d'ordre quatre H définissant l'écoulement plastique anisotrope défini par les six paramètres classiques de Hill F, G, H, L, M et N, et enfin les paramètres caractérisant l'évolution de l'endommagement ductile A, Y₀, γ et β.

Le multiplicateur plastique $\delta$ est déterminé par la condition de consistance associé au critère de chargédécharge ($f \le 0$; $\dot{\delta} \le 0$ et $\dot{\delta}f = 0$). Dans le domaine plastique, et après quelques développements algébriques le multiplicateur plastique est déterminé [1, 2] par :

$$ \dot{\delta} = \begin{cases} \frac{1}{h} \langle 3\mu\sqrt{1-D\tilde{n}:D} \rangle & \text{si } f = \dot{f} = 0 \\ 0 & \text{autre} \end{cases} \qquad (10) $$

où h est le module d'élastoplasticité strictement positif défini par :

$$ h = (2\mu\sqrt{1-D} + C)\tilde{n} : \tilde{n} = \tilde{n} - Q - a\tilde{n} : \frac{X}{\sqrt{1-D}} - b\frac{R}{\sqrt{1-D}} + \frac{\sigma_y}{2(1-D)^{\beta+1}}\left(\frac{Y-Y_0}{A}\right)^\gamma \quad (11) $$

## 2.2 Schéma de résolution globale

A partir du principe des puissances virtuelles, on peut définir par la méthode des éléments finis les équations discrétisées régissant le problème d'équilibre par :

$$ [M]\{\ddot{U}\} + \{F_{\text{int}}\} - \{F_{\text{ext}}\} = 0 \qquad (12) $$

où [M] est la matrice de masse cohérente de chaque élément (e), $\{F_{\text{int}}\}$ le vecteur forces intérieures, et $\{F_{\text{ext}}\}$ le vecteur forces extérieures incluant les forces de volume, de contact et de surface.

Le système algébrique (12), fortement non linéaire, exprime l'équilibre mécanique du travail à chaque pas de temps. Ce système est résolu en utilisant le solveur Dynamique Explicite (DE) d'Abaqus®. A chaque incrément de charge de la taille $\Delta t = t_{n+1} - t_n$ le procédé explicite d'intégration détermine les champs de déplacements suivants :

$$ \left\{ \begin{aligned} \{\ddot{U}\}_n &= [M]_n^{-1} \{\{F_{\text{ext}}\}_n - \{F_{\text{int}}\}_n\} \\ \{\dot{U}\}_{n+1/2} &= \{\dot{U}\}_{n-1/2} + \frac{\Delta t_{n+1} + \Delta t_n}{2} \{\ddot{U}\}_n \\ \{U\}_{n+1} &= \{U\}_n + \Delta t_{n+1} \{\dot{U}\}_{n+1/2} \end{aligned} \right. \qquad (13) $$

Les indices (n), (n+1/2) et (n+1) se rapportent au début, à la moitié et à la fin de l'incrément de temps. Notons que la taille de l'incrément de temps assure l'exactitude et la stabilité de la solution [4].

Le calcul des forces internes et des forces de contact exige le calcul à l'instant $t_{n+1}$ des contraintes de Cauchy $\underline{\sigma}_{n+1}$ et des efforts de contact.

## 2.3 Détermination des contraintes

Le calcul du tenseur de contraintes à l'instant $t_{n+1}$ nécessite l'intégration numérique des équations constitutives entièrement couplées du Tableau 1. Pour ce faire, l'Algorithme de Retour Radial, largement discuté dans la littérature [1, 2, 5], est appliqué.

La prédiction élastique, définie par l'indice « essai », est déterminée comme suit :

$$ \underline{\sigma}_{n+1}^{\text{essai}} = \underline{\sigma}_n + (1-D_n)[2\mu\Delta\underline{\epsilon} + \lambda\text{tr}(\Delta\underline{\epsilon})] \text{ sachant que } \Delta\underline{\epsilon} = \underline{D}\Delta t \text{ and } \underline{D} = \underline{L}^{\text{sym}} = (\underline{F} \ \underline{F}^{-1})^{\text{sym}} \qquad (14) $$

où $\underline{L}$ est le gradient de vitesse de déformation, $\underline{F} = 1 - \text{grad}(\Delta u)$ le gradient de transformation et $\Delta u$ l'incrément de déplacement appliqué à la structure durant l'intervalle de temps $[t_n, t_{n+1}]$. En appliquant la prédiction élastique (14) sur la variable tensorielle $\underline{Z} = \underline{\sigma} - \underline{X}$, nous obtenons la prédiction de fonction de charge définit ci dessous :
$$f^{\text{essai}} = \frac{1}{\sqrt{1-D_n}} (\|Z_n\| - R_n) - \sigma_y \quad (15)$$

Si $f^{\text{essai}} \le 0$ le comportement est purement élastique, les variables internes sont donc :

$$(\underline{\sigma}_{n+1} = \underline{\sigma}_n, \underline{X}_{n+1} = \underline{X}_n, R_{n+1} = R_n, \underline{\varepsilon}_{n+1}^p = \underline{\varepsilon}_n^p \text{ et } D_{n+1} = D_n).$$

Si $f^{\text{essai}} > 0$ le comportement est élastoplastique et une correction plastique est nécessaire sur les variables internes $\underline{\sigma}$, $\underline{X}$, $R$, $D$.

L'objectif de l'intégration numérique est alors de déterminer les inconnues $(\underline{\sigma}_{n+1}, \underline{X}_{n+1}, R_{n+1}, \underline{\varepsilon}_{n+1}^p, D_{n+1}$ et $\Delta\delta)$ à l'instant $t_{n+1}$ en appliquant la méthode de retour radial.

Une opération supplémentaire de réduction du nombre d'équations régissant le comportement va nous permettre un gain considérable en temps CPU [1, 2, 3, 5]. Après quelques manipulations algébriques [2, 3, 6] et en utilisant un algorithme asymptotique, on obtient ainsi le système réduit suivant :

$$ \begin{cases} f(\Delta\delta, \tilde{n}_{n+1}) = \frac{1}{\sqrt{1-D_{n+1}}} (\|Z_{n+1}\| - R_{n+1}) - \sigma_y = 0 \\ h(\Delta\delta, \tilde{n}_{n+1}) = \underline{H} : Z_{n+1} - \|Z_{n+1}\| \tilde{n}_{n+1} = 0 \end{cases} \quad \text{avec } \|Z_{n+1}\|_{\underline{H}} = \sqrt{Z_{n+1} : \underline{H} : Z_{n+1}} \quad (16) $$

$$ Z_{n+1} = Z_{n+1}^{\text{essai}} - \sqrt{1-D_{n+1}} \left[ 2\mu\Delta\delta + \frac{2}{3} \frac{C}{a} (1-e^{-a\Delta\delta}) \right] \tilde{n}_{n+1} \quad \text{ou} \quad Z_{n+1}^{\text{essai}} = (1-D_{n+1}) \left[ \frac{\sigma_{n+1}^{\text{essai}}}{1-D_n} - \frac{2}{3} C \underline{\alpha}_n e^{-a\Delta\delta} \right] \tilde{n}_{n+1} \quad (17) $$

La première équation (scalaire) du système (16) représente le critère de plasticité et la seconde (tensorielle) représente la normale $\tilde{n}$ à la fonction de charge. Ce système d'équations non-linéaires est résolu par la méthode itérative de Newton-Raphson.

## 2.4 Détermination des forces de contact

En mise en forme le contact/frottement joue un rôle important, à l'interface entre les outils et la pièce déformable. En conséquence, les forces de contact (avec ou sans le frottement) doivent être présentes dans les équations du mouvement. Cette prise en compte s'effectue à l'aide de la méthode de pénalité, qui est la plus répandue dans les codes de calculs par éléments finis, et qui consiste à rajouter un terme de pénalisation à l'énergie potentielle écrite sous la forme faible [4, 7, 8]. Dans le code Abaqus®, les conditions d'interface de contact sont prises en considération en utilisant l'approche des surfaces maître/esclave [4]. L'équation du frottement est basée sur le modèle de Coulomb, caractérisé par le coefficient de frottement $\mu$.

# 3 APPLICATION AU PROCEDE DE REFENDAGE EN 3D

## 3.1 Présentation du procédé

La méthodologie présentée ci-dessus est maintenant employée pour simuler numériquement le démarrage du refundage de tôles, avant enroulement sur la bobineuse, sur une longueur de 120mm. Les caractéristiques géométriques de la tôle sont donc 120mm*90mm*1mm et celles des outils D=80mm, largeur=10mm, plus petits qu'en réalité. Les entretoises (simples ou revêtues) ne sont pas modélisées.

<table><caption>Tableau 2. Paramètres du matériau utilisé</caption><tr><td rowspan="2">Elasticité</td><td>E</td><td>v</td><td colspan="3">σ<sub>y</sub></td></tr><tr><td>220 GPa</td><td>0.3</td><td colspan="3">410 MPa</td></tr><tr><td rowspan="4">Plasticité</td><td>a</td><td>b</td><td>C</td><td colspan="2">Q</td></tr><tr><td>31.5</td><td>6.5</td><td colspan="2">5300 MPa</td><td>2000 MPa</td></tr><tr><td>F</td><td>G</td><td>H</td><td colspan="2">L=M</td></tr><tr><td>0.43</td><td>0.46</td><td>0.53</td><td>1.5</td><td>1.63</td></tr><tr><td rowspan="2">Endommagement</td><td>Y<sub>0</sub></td><td>γ</td><td colspan="2">A</td><td>β</td></tr><tr><td>11</td><td>4.5</td><td colspan="2">28 MPa</td><td>2</td></tr></table>
Les simulations sont effectuées avec le code Abaqus/Explicit® en utilisant la formulation lagrangienne actualisée, le modèle de comportement mécanique étant implanté dans une routine VUMAT spécifique.

Le matériau de la tôle est supposé élastoplastique anisotrope endommageable, ses propriétés sont obtenues
à partir d'essais de cisaillement simple (voir Tableau 2).

La tôle est discrétisée avec 69800 éléments volumiques (éléments C3D8R à 8 nœuds avec intégration réduite). Une zone de largeur 3 mm a été choisie comme lieu de raffinement préalable du maillage (voir Figure 2). Les outils supposés indéformables sont discrétisés avec des éléments rigides (R3D4).

Les conditions aux limites ainsi que le chargement du problème sont appliqués en 2 étapes (Figure 2):

ETAPE 1 : Positionnement des couteaux

* Tôle : symétrie sur le plan (0,Y,Z),

* Outils : vitesse verticale de déplacement sur l'axe (OZ) $V_1$= 200mm/s sur les couteaux inférieurs et $V_1$= -200mm/s sur les couteaux supérieurs, jusqu'à une pénétration de 0,28mm de chaque couteau.

ETAPE 2 : Refendage proprement dit

* Tôle : symétrie sur le plan (0,X,Z) + vitesse de déplacement suivant (OY) de $V_2$=1000 mm/s,

* Outils : vitesse de rotation $W_1$=25 rd/s sur les couteaux inférieurs et $W_2$=-25rd/s sur les couteaux supérieurs, la rotation s'effectuant autour de l'axe des arbres cylindriques.

<figure><img src="image_2.png" /><figcaption>FIG. 1. Vue en sortie de refendeuse</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Maillage et conditions aux limites.</figcaption></figure>

**3.2 Résultats.**

La Figure 3 montre l'efficacité du modèle à reproduire la propagation de la fissure pendant le refendage. Concernant la bande découpée, on visualise sur cette même figure, d'une part la diminution des contraintes due à la redistribution mécanique au fur et à mesure de la découpe, et d'autre part son corollaire, à savoir l'apparition d'une déformation résiduelle. Cette déformation évolue pour se stabiliser lorqu'on atteint l'état stationnaire (pas tout à fait atteint ici). La figure 4 montre l'effort de coupe : partant d'une valeur non-nulle au point P1 (les couteaux ont pénétré dans la tôle), cet effort prend comme attendu une valeur constante (de l'ordre 18 kN) entre les point P2 et P3, puis décroît et s'annule en P4, lorsque la tôle est totalement cisaillée.

<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /><figcaption>FIG. 3. Répartition de la contrainte équivalente de Von Mises en fonction du déplacement des couteaux</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIG. 4. Evolution de l'effort de découpe en fonction du déplacement des couteaux</figcaption></figure>

## 4 Conclusion.

Dans cet article, nous avons illustré sur le procédé de refendage une méthodologie efficace pour modéliser numériquement les procédé industriels mettant en jeu de grandes déformations. Nous avons présenté d'une part les aspects théoriques et numériques du modèle, et d'autre part les résultats obtenus en ce qui concerne le refendage de tôles. On notera la capacité du modèle à réaliser sans difficulté la découpe, grâce à la prévision correcte de la localisation dans la tôle de l'endommagement induit par les couteaux.

Concernant la simulation du refendage, il s'agit ensuite d'intégrer l'ensemble des outils, et d'effectuer les calculs sur une longueur plus grande, permettant d'atteindre de façon sûre l'état stationnaire du procédé.

Au niveau du modèle, les perspectives envisagées au LASMIS sont les suivantes :

- introduction d'une formulation non-locale,

- prise en compte de l'endommagement des outils,

- prise en compte de l'anisotropie de l'endommagement.

## References

[1] Belamri, N (2004). Modélisation numérique du découpage de tôles laminées. Application au cisaillage et au refendage en 3D. Thèse de doctorat 3ème Cycle. Université de Technologie de Troyes.

[2] Saanouni K. & Chaboche J.L. (2003). Computational damage mechanics. Application to metal forming, Numerical and Computational methods, Ed. I. Miline, R.O. Ritchie and B. Karihaloo, Elsevier Oxford, Chap 7, vol. 3, pp. 321-376.

[3] Belamri N., Kkelifa K., Saanouni K., Cherouat A. and Autesserre P. (2004). Sheet metal blanking simulation using a coupled anisotropic plastic flow damage modelling. 7th International Conference ESAFORM, Trondheim, Norvège, Avril 28-30.

[4] Abaqus. (2004). Theory Manual. ABAQUS Version 6.4, HKS Inc.

[5] Simo J.C. and Hughes T.J.R. (1998). Computational Inelasticity. Springer-Verlag.

[6] Chaboche J-L., Cailletaud G. (1996). Integration methods for complex plastic constitutive equation. Computer Methods in Applied Mechanics and Engineering. vol.133, pp. 125-155.

[7] Laursen T.A. (2002). Computational contact mechanics. Springer-Verlag, Berlin

[8] Wriggers P and Panagiotopoulos P.D. (1999). New developments in contact problems. Springer-Verlag, Berlin Heidelberg.