# Développement d'une approche "Maximum Entropy" pour la modélisation thermo-mécanique du procédé de soudage par friction rotative

Mathieu Foca, Guillaume Racineux, Laurent Stainier

► To cite this version:

Mathieu Foca, Guillaume Racineux, Laurent Stainier. Développement d'une approche "Maximum Entropy" pour la modélisation thermo-mécanique du procédé de soudage par friction rotative. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722095

HAL Id: hal-01722095

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

Submitted on 2 Mar 2018

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 CC0 1.0 - Universal - International License
# Développement d'une approche "Maximum Entropy" pour la modélisation thermo-mécanique du procédé de soudage par friction rotative

Mathieu FOCA, Guillaume RACINEUX, Laurent STAINIER

GeM (UMR 6183 CNRS), Ecole Centrale de Nantes, {mathieu.foca, guillaume.racineux, laurent.stainier} @ec-nantes.fr

**Résumé** — Une modélisation multi-physique du procédé de soudage par friction rotative basée sur une approche Maximum Entropy (MaxEnt) est proposée. Cette approche permettra de résoudre des problèmes thermo-mécaniques couplés. Les déformations étant très importantes localement au voisinage de la soudure, le temps de remaillage lors d'une simulation par éléments finis classique est important. L'utilisation de cette méthode sans maillage à ce jour peu exploitée devrait permettre de réduire les temps de simulation ainsi que le phénomène de diffusion numérique.

**Mots clés** — Couplage thermo-mécanique, Maximum Entropy, Meshfree, Soudage par friction rotative

## 1 Introduction

Le soudage par friction rotative (RFW - Rotary Friction Welding) permet de réaliser des assemblages en bout de deux pièces dont l'une au moins a une symétrie de révolution. Il met en oeuvre l'effet thermique engendré dans leur plan de joint par la rotation rapide sous pression des pièces l'une contre l'autre. La puissance de chauffe résulte du couple résistant et de la vitesse de rotation. Du fait du fluage du matériau pendant la friction et lors du forgeage, la température de fusion n'est jamais atteinte. Ceci confère à la soudure des qualités exceptionnelles comparées à d'autres moyens de liaison tels que le soudage par résistance ou par étincelage, le soudage par faisceau d'électrons ou bien le soudage laser. Deux techniques de soudage par friction rotative sont principalement utilisées aujourd'hui : le soudage par friction pilotée, et le soudage par friction inertielle.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Principe du procédé de soudage par friction rotative, ici avec un cylindre plein.</figcaption></figure>

## 2 Modélisation

Pour obtenir une modélisation numérique cohérente du procédé, il convient de prendre en compte tous les aspects thermiques et mécaniques mis en jeu. Il s'agit donc de prendre en compte le couplage thermo-mécanique qui intervient lors de la soudure. En effet, les déformations sont très importantes et les températures atteignent des valeurs très élevées au niveau de la soudure. Il est donc important d'avoir un code de calcul qui puisse prendre en compte tout ces paramètres. L'ensemble de l'étude se fera en 2D
axisymétrique étant donnée la nature cylindrique (plein ou creux) des pièces à souder.

Afin d'évaluer la pertinence de l'étude et son efficacité, les résultats obtenus seront comparés aux résultats d'une simulation par éléments finis équivalente sur Abaqus, qui sera considérée comme étant la solution industrielle de référence. Pour une telle simulation, Abaqus fonctionne par remaillages successifs dès que la consommation matière atteint un cap prédéfini : par exemple, on ne fait pas plus de 0.5mm de consommation matière avec un seul maillage. Ainsi, on obtient une évolution du maillage telle que le montre la figure 2. Les deux inconvénients majeurs de cette méthode sont le temps de consacré au remaillage et la diffusion numérique qu'elle entraîne. En effet, ce deuxième aspect est important dans le contexte du travail présenté ici car l'objectif par la suite sera d'être capable de mener des étude métallurgique du matériau après soudage, ce qui demande une grande précision sur le champ thermique lors de la simulation.

La méthode Maximum Entropy a été préférée à d'autres méthodes telles que Smoothed Particuled Hydrodynamics (SPH) ou Element Free Galerkin car nous avons besoin ici de gérer des contacts et donc de parfaitement définir les contours géométriques de notre pièce à chaque instant.

L'objectif à plus long terme de ce travail est d'être capable de résoudre des problème de thermo-mécanique au sens large, des références bibliographiques [1] [2] montrant que la méthode donne de bons résultats sur des problèmes de mécanique ayant un comportement thermique adiabatique.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Évolution d'un maillage avec le code Abaqus [3]</figcaption></figure>

## 3 Aspect Expérimental

En parallèle du travail de simulation, des essais sur machine RDS60 sont effectués chez ACB afin de pouvoir par la suite vérifier nos résultats. Dans un premier temps, on soude des pièces en acier P295GH. Ces pièces ont les géométries suivantes : tube épais, tube mince et cylindre plein. Quand l'influence de tous les paramètres de soudage sera bien maitrisée, des pièces en alliage de titane TA6V seront utilisées.

Le procédé de soudage par friction rotative est divisé en quatre phases :

- Phase 1 : Rampe de montée en force

- Phase 2 : Friction à effort et vitesse de rotation constants.

- Phase 3 : Transitoire vers effort de forge imposé.

- Phase 4 : Refroidissement sous maintien en force.

Les paramètres d'entrée de la machine sont les suivants, les valeurs numériques indiquées étant celles correspondant à notre premier plan d'essai :
- La force exercée sur la pièce en rotation. Celle-ci varie entre 35 kN et 155 kN, soit entre 6% et 26% de la capacité de fonctionnement de la machine.

- La vitesse de rotation V : elle varie entre 600 tr/min et 900 tr/min, soit entre 66% et 100% de la capacité de fonctionnement de la machine.

- La consommation matière U : quantité de matière qui sera éjectée tout au long de l'opération et qui constituera le bourrelet au niveau de la zone de soudure. Elle varie entre 4 mm et 8 mm.

L'objectif est de comparer les aspects suivant entre simulation et expérimental :

- Temps d'échauffement

- Vitesse de consommation matière

- Champ de température

- Effet de la pression de forge

- Taille de la Zone Affectée Thermiquement (ZAT)

L'utilisation d'une caméra thermique permet en outre d'obtenir la température de surface au voisinage de la soudure.

# 4 La méthode Maximum Entropy

## 4.1 Origines

La méthode MaxEnt utilise la notion d'entropie au sens de la théorie de l'information comme la définit Shannon [4]. Dans cette théorie, l'entropie est définie comme une grandeur permettant de quantifier l'incertitude. Chercher une solution maximisant l'entropie revient à chercher le cas où le résultat est le plus incertain. Illustrons cet aspect avec un exemple. Soit un ensemble d'événements A = {A₁, ..., Aₙ} associés aux probabilités {p₁, ..., pₙ}. Considérons les deux cas suivant :

$$B_1 = \begin{pmatrix} A_1 & A_2 \\ 0.5 & 0.5 \end{pmatrix} \text{ et } B_2 = \begin{pmatrix} A_1 & A_2 \\ 0.9 & 0.1 \end{pmatrix}$$

L'issue est évidemment plus incertaine dans le cas B₁ que dans le cas B₂ (où le résultat sera très certainement A₁). L'entropie est donc plus élevée dans le cas B₁. En ce basant sur trois critères mathématiques, Shannon [4] établit la formulation de l'entropie H par la formule suivante :

$$H(B) = H(p_1, \dots, p_n) = - \sum_{a=1}^{n} p_a \log p_a \quad (1)$$

On a donc :

$$H(B_1) = 0.301 \text{ et } H(B_2) = 0.141$$

## 4.2 Mise en place du problème

Considérons un ensemble de nœuds distincts X = {xₐ, a = 1, ..., N} ⊂ ℝ<sup>d</sup>. L'espace convexe généré par X est donné par la formule suivante [5] :

$$\mathrm{conv}X = \{ \mathbf{x} \in \mathbb{R}^d | \mathbf{x} = \mathbf{X}\mathbf{\lambda}, \mathbf{\lambda} \in \mathbb{R}_+^N, \mathbf{1} \cdot \mathbf{\lambda} = 1 \} \quad (2)$$

où ℝ<sup>N</sup> est l'orthant positif, 1 ∈ ℝ<sup>N</sup> est le vecteur rempli de 1 et X est une matrice d x N dont les colonnes sont les coordonnées des vecteurs positions des nœuds de X. X étant fini, alors convX est convexe et compact.

Soit u : convX → ℝ dont les évaluations {uₐ = u(xₐ), a = 1, ..., N} sont connues sur X. On veut construire des approximations de u (qui peut par exemple représenter un champ de déplacement) de la forme suivante :

$$u^h(\mathbf{x}) = \sum_{a=1}^{n} p_a(\mathbf{x})u_a \quad (3)$$
où les fonctions $p_a : \text{conv}X \to \mathbb{R}$ sont les fonctions de forme ou fonctions d'interpolation. Ces fonctions doivent répondre aux conditions suivantes :

$$ \sum_{a=1}^{n} p_a(\mathbf{x}) = 1 \qquad (4) $$

$$ \sum_{a=1}^{n} p_a(\mathbf{x})\mathbf{x}_a = \mathbf{x} \qquad (5) $$

A partir de ces relations et de la notion d'entropie établie par Shannon, Arroyo et Ortiz [1] montrent qu'il est possible de faire le lien entre les deux. En supposant que le choix des fonctions de forme doit être le moins biaisé possible et d'après le principe d'entropie maximale de Jaynes [6], il convient de résoudre le problème *Maximum Entropy* (ME) suivant :

$$ \begin{array}{ll} \text{(ME) Maximiser} & H(\mathbf{p}) = - \sum_{a=1}^{n} p_a \log p_a \\ \text{tel que} & p_a \ge 0, a = 1, \dots, N \\ & \sum_{a=1}^{n} p_a(\mathbf{x}) = 1 \\ & \sum_{a=1}^{n} p_a(\mathbf{x})\mathbf{x}_a = \mathbf{x} \end{array} $$

De plus, afin de contrôler le degré de corrélation entre les fonction d'interpolation en $x$ et les valeurs aux nœuds proches, il convient d'ajouter une notion de localité, ce qui se traduit par le problème (RAJ) [7] :

$$ \begin{array}{ll} \text{(RAJ) À x fixé, minimiser} & U(\mathbf{x}, \mathbf{p}) = \sum_{a=1}^{n} p_a |\mathbf{x} - \mathbf{x}_a|^2 \\ \text{tel que} & p_a \ge 0, a = 1, \dots, N \\ & \sum_{a=1}^{n} p_a(\mathbf{x}) = 1 \\ & \sum_{a=1}^{n} p_a(\mathbf{x})\mathbf{x}_a = \mathbf{x} \end{array} $$

Pour obtenir le compromis optimal entre ces deux problèmes, il faut donc résoudre le problème suivant :

$$ \begin{array}{ll} \text{(LME)}_{\beta} \text{ À x fixé, minimiser} & f_{\beta}(\mathbf{x}, \mathbf{p}) = \beta U(\mathbf{x}, \mathbf{p}) - H(\mathbf{p}) \\ \text{tel que} & p_a \ge 0, a = 1, \dots, N \\ & \sum_{a=1}^{n} p_a(\mathbf{x}) = 1 \\ & \sum_{a=1}^{n} p_a(\mathbf{x})\mathbf{x}_a = \mathbf{x} \end{array} $$

Dans leur article, Arroyo et Ortiz [1] montrent également que la solution au problème $(\text{LME})_{\beta}$ est unique et s'écrit :

$$ p_{\beta a} = \frac{1}{Z(\mathbf{x}, \mathbf{\lambda}^*(\mathbf{x}))} \exp[-\beta|\mathbf{x} - \mathbf{x}_a|^2 + \mathbf{\lambda}^*.(\mathbf{x} - \mathbf{x}_a)] \qquad (6) $$

où :

$$ Z(\mathbf{x}, \mathbf{\lambda}) = \sum_{a=1}^{n} \exp[-\beta|\mathbf{x} - \mathbf{x}_a|^2 + \mathbf{\lambda}.(\mathbf{x} - \mathbf{x}_a)] \qquad (7) $$

$$ \mathbf{\lambda}^*(\mathbf{x}) = \arg \min_{\mathbf{\lambda} \in \mathbb{R}^d} \log Z(\mathbf{x}, \mathbf{\lambda}) \qquad (8) $$

Le paramètre $\lambda$ correspond aux multiplicateurs de Lagrange imposant les conditions (4) et (5).

Le paramètre $\beta$ est à fixer au cours de la simulation en fonction de la richesse de l'interpolation qu'on souhaite obtenir. $\beta = 0$ signifie que le support des fonctions de forme sera le "maillage" en entier alors que $\beta = +\infty$ signifie qu'on néglige l'aspect biaisé du choix des fonctions de forme. On montre que quelque soit $\beta \in [0, +\infty[$ et quelque soit $\mathbf{x} \in \text{conv}X$, la solution existe et est unique. En pratique, on utilise également $\gamma$, qui est lié à $\beta$ par la relation :

$$ \gamma = h^2 \beta \qquad (9) $$

où $h$ est la longueur caractéristique entre deux nœuds.
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Fonctions de forme en 1D pour différentes valeurs de $\gamma$, $\gamma = \{0, 0.5, 1.8, 5\}$.</figcaption></figure>

# 5 Travaux et premiers résultats

Dans la bibliographie, la méthode MaxEnt n'a été utilisée que pour résoudre des problèmes de mécanique [1]. Avant de résoudre des problèmes thermo-mécanique fortement couplés, on vérifie que la méthode donne de bons résultats sur des problèmes thermiques simples. Pour ce faire, les routines Matlab partagées par Arroyo sont utilisées comme base.

## 5.1 Visualisation des fonctions de forme

### 5.1.1 Fonctions de forme en 1D

Dans cet exemple 1D, les nœuds sont au nombre de cinq et ont les coordonnées suivantes :

$$x_a = 0.01a, \quad \forall a \in [0; 1; 2; 3; 4]$$

Dans ce cas, $h = 0.01$. Chacune des cinq fonctions est évaluée en chaque nœud. On donne à la figure 3 leurs représentations pour différentes valeurs de $\gamma$.

Les valeurs données dans le tableau 5.1.1 permettent de vérifier qu'on respecte bien les conditions (4) et (5).

### 5.1.2 Fonctions de forme en 2D

De la même façon, on représente sur la figure 4 les fonctions de forme sur un modèle 2D.
<table><caption>Tableau 1 – Valeurs des fonctions de forme sur un exemple 1D, $\gamma = 1.8$.</caption><thead><tr><th></th><th>Nœud 1</th><th>Nœud 2</th><th>Nœud 3</th><th>Nœud 4</th><th>Nœud 5</th></tr></thead><tbody><tr><th>Fonction 1</th><td>1.000</td><td>1.247e<sup>-1</sup></td><td>5.605e<sup>-4</sup></td><td>6.828e<sup>-8</sup></td><td>1.562e<sup>-27</sup></td></tr><tr><th>Fonction 2</th><td>4.401e<sup>-5</sup></td><td>7.511e<sup>-1</sup></td><td>1.241e<sup>-1</sup></td><td>5.558e<sup>-4</sup></td><td>1.740e<sup>-18</sup></td></tr><tr><th>Fonction 3</th><td>5.249e<sup>-11</sup></td><td>1.236e<sup>-1</sup></td><td>7.507e<sup>-1</sup></td><td>1.236e<sup>-1</sup></td><td>5.294e<sup>-11</sup></td></tr><tr><th>Fonction 4</th><td>1.740e<sup>-18</sup></td><td>5.558e<sup>-4</sup></td><td>1.241e<sup>-1</sup></td><td>7.511e<sup>-1</sup></td><td>4.401e<sup>-5</sup></td></tr><tr><th>Fonction 5</th><td>1.562e<sup>-27</sup></td><td>6.828e<sup>-8</sup></td><td>5.605e<sup>-4</sup></td><td>1.247e<sup>-1</sup></td><td>1.000</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Fonctions de formes sur un modèle 2D</figcaption></figure>

## 5.2 Résolution de problèmes thermiques simples

La résolution d'un problème thermique passe par la résolution du système suivant :

$$[C]\{\dot{T}\} + [K]\{T\} = \{q_e\} \qquad (10)$$

avec $T$ le vecteur contenant les températures nodales, $C$ la matrice de capacité thermique, $K$ la matrice de conductibilité et $q_e$ le vecteur contenant les flux nodaux imposés.

### 5.2.1 Calcul de la matrice de conductibilité $K$

Le calcul de la matrice $K$ est donné par l'équation suivante :

$$\forall (a, b), K_{a,b} = \sum_{a=1}^{n} k \nabla p_{\beta a} \cdot \nabla p_{\beta b} dV \qquad (11)$$

avec $k$ la conductivité thermique du matériau, $\nabla p_{\beta i}$ le gradient de la fonction de forme au nœud $i$, obtenu en dérivant la fonction de forme $p_{\beta i}$ et $\Omega_e$ le convexe formé par l'ensemble des nœuds.

### 5.2.2 Calcul de la matrice de capacité thermique $C$

Le calcul de la matrice $C$ est donné par l'équation suivante :

$$\forall (a, b), C_{ab} = \sum_{a=1}^{n} \rho c_v p_{\beta a} p_{\beta b} dV \qquad (12)$$
avec $\rho$ la masse volumique du matériau, la chaleur spécifique du matériau, $p_{\beta i}$ la valeur de la fonction de forme au nœud i et $\Omega_e$ le convexe formé par l'ensemble des nœuds.

Afin de pouvoir comparer avec les résultats fournis par Abaqus, il convient d'utiliser la même matrice C. Or, Abaqus utilise une matrice C' condensée sur la diagonale donnée la formule suivante. On utilisera donc la version condensée de cette matrice pour notre calcul.

$$ \forall (i, j) \in [1, n]^2, C'_{ij} = \delta_{ij} \sum_{k=1}^{n} C_{ik} \qquad (13) $$

### 5.2.3 Schéma d'intégration numérique

On résout le problème thermique en intégrant dans le temps par différence finie. En particulier, la méthode des trapèzes généralisés est utilisée :

$$ \begin{cases} [C']\{\dot{T}_{n+1}\} + [K]\{T_{n+1}\} = \{q_e(t_{n+1})\} \\ \{T_{n+1}\} = \{T_n\} - (1-\alpha)\Delta t\{\dot{T}_n\} + \alpha\Delta t\{\dot{T}_{n+1}\} \end{cases} \qquad (14) $$

avec :

* $\alpha = 0$ : schéma explicite

* $\alpha = 1/2$ : schéma de Crank-Nicholson

* $\alpha = 1$ : schéma Euler-implicite

Pour la suite, le schéma de Crank-Nicholson est choisi.

## 5.3 Cas test 1D, comparaison MaxEnt/FEM

Afin de vérifier la validité de notre modèle, on compare sur des problèmes très simple les résultats obtenus avec les deux méthodes. On ne présentera ici qu'un des cas tests réalisés.

Soit une barre composée d'un matériau homogène de longueur $L = 0.03m$. Le maillage éléments finis et la position des nœuds pour MaxEnt sont donnés par la figure 5. Dans cet exemple, on notera $T_0 = 273K$. Initialement, on fixe la température de la plaque à $T = T_0$ puis $\forall t \ge 0$, on impose la température sur un bord égale à $T = 400K$, ce qui se traduit mathématiquement ainsi : Conditions initiales :

$$ T(x, 0) = 0, \forall x \in [0, L] $$

$$ \dot{T}(x, 0) = 0, \forall x \in [0, L] $$

Conditions limites :

$$ T(0,t) = T_0, \forall t \ge 0 $$

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Maillage MEF (l) et nœuds et points matériels MaxEnt (●) pour l’exemple de plaque 1D</figcaption></figure>

La figure 6 compare les résultats obtenus avec une méthode élément classique avec Abaqus (traits pleins) et avec la méthode MaxEnt (traits discontinus). Les résultats sont très similaires, on valide donc la méthode pour la résolution de problèmes thermiques.

## 6 Conclusion

Dans ce travail, on propose d'appliquer le principe d'entropie maximale à la résolution de problème thermo-mécanique fortement couplés. La cohérence sur des problèmes mécaniques [1] et thermiques simples nous encourage sur cette voie. Après implémentation du code permettant la simulation numérique de tels procédés, nous devrions être en mesure de l'appliquer à la modélisation du procédé de soudage par friction rotative. Celle-ci sera suffisamment fine pour pouvoir par la suite faire des études métallurgiques afin de déterminer à l'avance les paramètres de soudage présentés dans ce papier.
<figure><img src="image_7.png" /><figcaption>Fig. 6 – Comparaison MEF (—) / MaxEnt (– · –) sur un problème de propagation thermique simple</figcaption></figure>

Remerciements

Ce travail a été réalisé dans le cadre du projet Friction Rotative à Entrainement Direct (FRED), financé par OSEO, en collaboration avec ACB, l'ENSAM (LAMPA) et Polytech' Nantes (IMN).

Références

[1] M. Arroyo, M. Ortiz. *Local maximum-entropy approximation schemes : a seamless bridge between finite elements and meshfree methods*, International journal for numerical methods in engineering, 65 :2167-2202, 2006.

[2] B. Li, A. Kidane, G. Ravichandran, M. Ortiz. *Verification and validation of the Optimal Transportation Meshfree (OTM) simulation of terminal ballistics*, INT J IMPACT ENG, 42 :25-36, 2012.

[3] Abaqus Example Problems Manual, 1.3.18 Inertia welding simulation using Abaqus/Standard and Abaqus/CAE, SIMULIA, Abaqus 6.10.

[4] C. E. Shannon. *A Mathematical Theory of Communication*, The Bell System Technical Journal, Éditeur, 27 :379-423, 27 :623-656, 1948.

[5] R. T. Rockafellar. *Convex Analysis*, Princeton University Press, Princeton NJ, 1970.

[6] E. T. Jaynes. *Information theory and statistical mechanics*, Physical Review, **106**(4):620-630, 1957.

[7] V. T. Rajan. *Optimality of the Delaunay triangulation*, Discrete and Computational Geometry, **12**(2):189-202, 1994.