# Modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température

Thomas Heuzé, J.B. Leblond, Jean-Michel Bergheau

► To cite this version:

Thomas Heuzé, J.B. Leblond, Jean-Michel Bergheau. Modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température. 20e Congrès Français de Mécanique, Aug 2011, Besançon, France. hal-02387437

HAL Id: hal-02387437

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

Submitted on 29 Nov 2019

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 CC BY 4.0 - Attribution - International License
# Modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température

T. Heuzéª, J.B. Leblondª, J.M. Bergheauᵇ

a. UPMC, Institut Jean Le Rond d'Alembert, 4 place jussieu tour 55-65, 75252 PARIS Cedex 05

b. ENISE, LTDS, 58 rue Jean Parot 42023 Saint Étienne Cedex 2

## Résumé :

On développe dans ces travaux un outil numérique permettant de simuler le procédé Friction Stir Spot Welding. Le modèle est basé sur un couplage fluide/solide permettant de décrire l'état de la matière dans toute la structure. Ce modèle s'appuie sur l'élément fini P1+/P1, développé dans le cas d'un couplage thermomécanique fort, à la fois en mécanique des fluides et en mécanique des solides. Le couplage est associé avec une approche euléro-lagrangienne et intégré au sein d'une nouvelle option du code SYSWELD®. On présente ici une première simulation du procédé Friction Stir Spot Welding.

## Abstract :

A numerical tool for the simulation of the Friction Stir Spot Welding process is developed. The model is based on a fluid/solid coupling which describes the state of the material in the structure. The model includes the P1+/P1 finite element, developed with a strong thermomechanical coupling, both for fluid mechanics and solid mechanics. The coupling, associated with an Arbitrary Eulerian Lagrangian approach, is implemented in a new option of the finite element code SYSWELD®. We present here a first simulation of the Friction Stir Spot Welding process.

**Mots clefs :** Élément P1+/P1 ; Friction Stir Spot Welding ; Couplage fluide/solide

## 1 Introduction

Le procédé de soudage ponctuel par friction et malaxage (Friction Stir Spot Welding [4]) est un procédé dérivé du soudage par friction et malaxage (Friction Stir Welding [6]) consistant à créer un joint de soudure « ponctuel » entre deux tôles superposées par pénétration d'un outil en rotation dans la matière. L'échauffement généré par le frottement ainsi que le mouvement de la matière entraînée génère deux phases, la première au voisinage de l'outil atteignant un état pâteux, l'autre restant solide.

Une approche eulérienne lagrangienne arbitraire (ALE), dans laquelle le maillage suit la matière dans la phase solide mais pas dans la zone pâteuse, permet de simuler de nombreux tours de l'outil dans la matière tout en suivant les bords des tôles. De plus, des modèles de comportement cohérents avec l'état de la matière dans chacune des zones peuvent être inclus, aussi la partie pâteuse est modélisée comme un fluide visqueux, les effets d'histoire étant pris en compte dans la partie solide.

Le modèle numérique développé est basé sur l'élément fini P1+/P1 [1]. Cet élément mixte est développé avec une formulation température/vitesse/pression en mécanique des fluides (pour un écoulement laminaire incompressible transitoire [2, 3]) et en mécanique des solides dans le cadre des grandes transformations.

Ces deux éléments sont couplés au sein d'une nouvelle option créée dans le code SYSWELD® [5], et présentés dans la section 2. La transition fluide/solide est effectuée au moyen d'un test explicite sur une température moyenne par élément, l'interface passe alors par les faces des éléments du maillage. Une procédure d'actualisation de la géométrie associée à l'approche euléro-lagrangienne est effectuée à convergence. Une première simulation simplifiée du procédé Friction Stir Spot Welding est présentée en section 3.
# 2 Couplage fluide/solide

## 2.1 Comportements thermomécaniques fluide et solide

La partie pâteuse est modélisée comme un fluide visqueux incompressible et décrite avec une loi de Norton-Hoff. Le solide est modélisé en grandes transformations, on postule une partition additive du taux de déformation eulérien **D** faisant donc intervenir une loi d'hypoélasticité. Les équations (1) résument les équations constitutives associées aux comportements mécaniques fluide et solide.

$$
\begin{array}{l@{\qquad}l}
\begin{array}{c}
\text{Fluide} \\
\begin{array}{r@{\,}l}
\boldsymbol{\sigma} & = p\mathbf{1} + \mathbf{s} \\
\mathbf{q} & = -k\nabla T
\end{array}
\end{array}
&
\begin{array}{c}
\text{Solide} \\
\begin{array}{r@{\,}l}
\mathbf{D} & = \mathbf{D}^e + \mathbf{D}^p + \varepsilon^{th}\mathbf{1} \\
\mathbf{s} & = \mathcal{F}(\mathbf{v}(\tau), \tau \le t) \\
\mathbf{D}^e & = \frac{1+\nu}{E}\nabla\boldsymbol{\sigma} - \frac{\nu}{E}(\text{tr}\,\nabla\boldsymbol{\sigma})\mathbf{1} \\
\text{tr } \mathbf{D}^p & = 0
\end{array}
\end{array}
\end{array}
\tag{1}
$$

où $\nabla\boldsymbol{\sigma}$ est une dérivée objective des contraintes de Cauchy, $p$ la pression hydrostatique, $\mathbf{s}$ le déviateur des contraintes, $\mathbf{q}$ le flux surfacique de chaleur, $\mu$ la viscosité dynamique du fluide, $E$ et $\nu$ le module d'Young et le coefficient de Poisson, $k$ la conductivité thermique, $\varepsilon^{th}$ le taux de déformation « thermométallurgique » provenant des dilatations thermiques et des changements de volume associés à la métallurgie, et $\mathcal{F}$ est une fonctionnelle décrivant la relation de comportement dépendant de l'histoire du solide, elle peut dépendre ou non de la vitesse de chargement. Les paramètres matériau dépendent à priori de la température.

## 2.2 Forme faible du problème

Le problème à résoudre fait intervenir les équations d'énergie, de mouvement et les liaisons internes de chaque zone. La formulation mixte du problème s'écrit sous la forme suivante, en notant $\chi$ les coordonnées associées à la configuration ALE :

$$
\begin{align}
& \text{Trouver } (T, \mathbf{v}, p) \in (\mathcal{T}_{ad} \times \mathcal{V}_{ad} \times \mathcal{P}_{ad}), \forall t \in [0, T], \text{ tel que } \forall (T^*, \mathbf{v}^*, p^*) \in (\mathcal{T}_{ad}^0 \times \mathcal{V}_{ad}^0 \times \mathcal{P}_{ad}), \nonumber \\
& \left\|
\begin{array}{@{}l@{}}
  - \int_{\Omega_{\chi}} k \nabla T \cdot \nabla T^* \, d\Omega + \int_{\partial\Omega_q} \phi^d T^* \, dS + \int_{\Omega_{\chi}} \left\{ \begin{array}{@{}l@{}}
    \boldsymbol{\sigma} : \mathbf{D} \\
    \beta \boldsymbol{\sigma} : \mathbf{D}^p
  \end{array} \right\} T^* \, d\Omega = \int_{\Omega_{\chi}} \rho C \left( \frac{\partial T}{\partial t} + \begin{pmatrix} \nabla T \cdot \mathbf{c} \\ 0 \end{pmatrix} \right) T^* \, d\Omega \\
  - \int_{\Omega_{\chi}} \left\{ \begin{array}{@{}l@{}}
    2\mu \mathbf{D} \\
    \mathbf{s}(\mathbf{v}(\tau), \tau \le t)
  \end{array} \right\} : \mathbf{D}^* \, d\Omega - \int_{\Omega} p \operatorname{div} \mathbf{v}^* \, d\Omega + \int_{\partial\Omega_F} \mathbf{F}_d \cdot \mathbf{v}^* \, dS = \int_{\Omega_{\chi}} \rho \left( \begin{pmatrix} \frac{\partial \mathbf{v}}{\partial t} + \nabla \mathbf{v} \cdot \mathbf{c} \\ \mathbf{0} \end{pmatrix} \right) \mathbf{v}^* \, d\Omega \tag{2} \\
  \int_{\Omega_{\chi}} p^* \left( \operatorname{div} \mathbf{v} - \begin{pmatrix} 0 \\ \frac{\dot{p}}{\kappa} + 3\varepsilon^{th} \end{pmatrix} \right) \, d\Omega = 0 \\
  T(\mathbf{x}, t = 0) = T_0 \\
  \mathbf{v}(\mathbf{x}, t = 0) = \mathbf{v}_0
\end{array}
\right. \nonumber
\end{align}
$$

où $\mathbf{c} = \mathbf{v} - \hat{\mathbf{v}}$ est la vitesse convective définie à partir de la vitesse matérielle $\mathbf{v}$ et de la vitesse du maillage $\hat{\mathbf{v}}$, $\phi^d$ et $\mathbf{F}_d$ représentent le flux surfacique de chaleur et les tractions surfaciques imposés, $\rho$ la masse volumique, $C$ la capacité calorifique, $\beta$ le coefficient de Quinney-Taylor généralement fixé à 0.9 et $\kappa$ le module de compressibilité élastique. Les accolades dans les équations permettent de distinguer les termes associés aux phases fluide (en haut) et solide (en bas), cette représentation symbolique sera conservée par la suite.

## 2.3 L'élément fini P1+/P1

Le tétraèdre P1+/P1 [1], servant de base au couplage, est un élément fini linéaire enrichi (Figure 1) assurant, dans la version proposée, la continuité des champs de température, vitesse et pression. Le
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Tétraèdre P1+/P1 de référence.</figcaption></figure>

champ bulle ajouté dans l'interpolation de la vitesse permet à cet élément de pouvoir traiter, au travers d'une formulation mixte, une liaison interne de type incompressibilité.

Le champ de vitesse est alors interpolé sur l'élément comme suit :

$$ \mathbf{v}^h(\mathbf{x}) = N_p(\xi, \eta, \zeta)\mathbf{v}^{(p)} + N_b(\xi, \eta, \zeta)\boldsymbol{\lambda} \quad 1 \le p \le 4 \quad (3) $$

où $p$ est le numéro d'un noeud sommet de l'élément, $b$ l'indice concernant le noeud bulle et les $\boldsymbol{\lambda}$ sont les degrés de libertés liés au noeud bulle, homogènes à des vitesses.

## 2.4 Formulation éléments finis

La discrétisation par éléments finis de la forme faible (2) conduit à l'écriture du système non linéaire d'équations semi-discrètes suivant :

$$ M\dot{\mathbf{q}} + \left\{ \begin{matrix} \mathbf{f}^{\text{conv}} \\ \mathbf{0} \end{matrix} \right\} + \mathbf{f}^{\text{int}} = \mathbf{f}^{\text{ext}} \quad (4) $$

le vecteur degré de liberté $\mathbf{q}$ étant défini tel que $\mathbf{q}^T = \{\mathbf{T} \ \mathbf{v} \ \mathbf{p} \ \boldsymbol{\lambda}\}$. Les efforts internes $\mathbf{f}^{\text{int}}$, externes $\mathbf{f}^{\text{ext}}$, et les efforts dûs aux termes de convection $\mathbf{f}^{\text{conv}}$ sont définis, en détaillant les composantes associées au $p$-ième noeud sommet ($p = 1$ à 4) et au noeud bulle ($b$), de la façon suivante :

$$ \begin{aligned} \mathbf{f}^{\text{int}} &= \sum_{e=1}^{N_{el}} \left( \begin{matrix} f_{T}^{\text{int},(p)} = - \int_{\Omega_e} \nabla N^{(p)} \cdot \mathbf{q} \, d\Omega \\ \mathbf{f}_{v}^{\text{int},(p)} = \int_{\Omega_e} \mathbf{B}_{p}^{T} \, \mathrm{sd}\Omega + \int_{\Omega_e} \nabla N^{(p)} \, p \, d\Omega \\ f_{p}^{\text{int},(p)} = \int_{\Omega_e} N^{(p)} \left\{ \begin{matrix} \text{div } \mathbf{v} \\ \text{div } \mathbf{v} - \frac{\dot{p}}{\kappa} - 3\dot{\varepsilon}^{\text{th}} \end{matrix} \right\} d\Omega \\ \mathbf{f}_{b}^{\text{int}} = \int_{\Omega_e} \mathbf{B}_{b}^{T} \, \mathrm{sd}\Omega + \int_{\Omega_e} \nabla N^{(b)} \, p \, d\Omega \end{matrix} \right) \\ \mathbf{f}^{\text{ext}} &= \sum_{e=1}^{N_{el}} \left( \begin{matrix} f_{T}^{\text{ext},(p)} = \int_{\partial\Omega_e \cap \partial\Omega_q} \phi^d N^{(p)} \, dS + \int_{\Omega_e} \left\{ \begin{matrix} \boldsymbol{\sigma} : \mathbf{D} \\ \beta \boldsymbol{\sigma} : \mathbf{D}^p \end{matrix} \right\} N^{(p)} \, d\Omega \\ \mathbf{f}_{v}^{\text{ext},(p)} = \int_{\partial\Omega_e \cap \partial\Omega_F} \mathbf{F}_d N^{(p)} \, dS \\ 0 \\ \mathbf{f}_{b}^{\text{ext}} = \mathbf{0} \end{matrix} \right) \\ \mathbf{f}^{\text{conv}} &= \sum_{e=1}^{N_{el}} \left( \begin{matrix} f_{T}^{\text{conv},(p)} = \int_{\Omega_e} \rho C(\nabla T) \cdot \mathbf{c} \, \bar{N}^{(p)} \, d\Omega \\ \mathbf{f}_{v}^{\text{conv},(p)} = \int_{\Omega_e} \rho(\nabla \mathbf{v}) \cdot \mathbf{c} \, N^{(p)} \, d\Omega \\ 0 \\ \mathbf{f}_{b}^{\text{conv}} = \int_{\Omega_e} \rho(\nabla \mathbf{v}) \cdot \mathbf{c} \, N^{(b)} \, d\Omega \end{matrix} \right) \end{aligned} \quad (5) $$
où $\mathbf{B}_p$ est la matrice contenant les gradients de la fonction de forme associée au noeud $p$, reliant les taux de déformation aux composantes de la vitesse du noeud $p$. On rappelle que $\int \mathbf{F}_d N^{(b)} dS = \mathbf{0}$ car la fonction bulle s'annule sur le pourtour de l'élément.

Un algorithme aux différences finies de type Euler implicite est utilisé pour la résolution. Dans le cas de la méthode de Newton-Raphson, la linéarisation des équations discrètes effectuée à l'itération $k$ du calcul conduit à la résolution du système matriciel suivant :

$$ \mathbf{M}^{(k)} \dot{\mathbf{q}}^{(k+1)} + \mathbf{K}^{(k)} \delta \mathbf{q}^{(k)} = \mathbf{R}^{(k)} \quad (6) $$

où l'incrément $\delta \mathbf{q}^{(k)}$ est défini par $\delta \mathbf{q}^{(k)} = \mathbf{q}^{(k+1)} - \mathbf{q}^{(k)}$. Les matrices de rigidité $\mathbf{K}^{(k)}$ et le vecteur résidu $\mathbf{R}^{(k)}$ sont définis comme suit :

$$ \mathbf{R}^{(k)} = \mathbf{f}^{\text{ext}} - \begin{Bmatrix} \mathbf{f}^{\text{conv}} \\ \mathbf{0} \end{Bmatrix} - \mathbf{f}^{\text{int}} \\ \mathbf{K}^{(k)} = - \frac{\partial \mathbf{R}}{\partial \mathbf{q}} \Big|^{(k)} \quad (7) $$

On notera que $\mathbf{R}^{(k)}$ ne s'annule pas à la convergence des itérations. Cette présentation est cohérente avec la structure du code : la routine élémentaire calcule les matrices de masse et de rigidité et le vecteur résidu, et la discrétisation temporelle des termes d'inertie est effectuée par une autre routine dans la boucle sur les mailles. Pour gagner du temps lors de la résolution, le noeud bulle est condensé dans les équations relatives aux noeuds sommets [3]. La résolution à chaque itération est alors effectuée sur l'incrément du vecteur des degrés de liberté $\bar{\mathbf{q}}$ réduit après condensation :

$$ \mathbf{A}^{(k)} \delta \bar{\mathbf{q}}^{(k)} = \mathbf{B}^{(k)} \quad \text{avec} \quad \left| \begin{array}{l} \mathbf{A}^{(k)} = \left( \frac{\bar{\mathbf{M}}^{(k)}}{\Delta t} + \bar{\mathbf{K}}^{(k)} \right) \\ \mathbf{B}^{(k)} = \bar{\mathbf{R}}^{(k)} - \bar{\mathbf{M}}^{(k)} \dot{\bar{\mathbf{q}}}^{(k)} \end{array} \right. \quad (8) $$

## 2.5 Actualisation de la géométrie

La procédure d'actualisation du maillage de la partie fluide est effectuée à la convergence du calcul, elle se déroule en trois temps : (i) les noeuds sont codés en fonction de la zone dans laquelle ils se trouvent (Figure 2(a)), (ii) la position des noeuds de la surface libre est actualisée à partir de la projection normale de la vitesse matérielle convergée (Figure 2(b)) et (iii) les noeuds intérieurs de la zone fluide sont replacés de la façon la plus régulière possible par une technique de replacement barycentrique (Figure 2(c)).

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Actualisation de la géométrie</figcaption></figure>

## 3 Simulation numérique du Friction Stir Spot Welding

On présente une simulation simplifiée du Friction Stir Spot Welding : (i) le contact n'est pas pris en compte aussi la simulation est menée sur une structure monobloc, (ii) le contact frottant est modélisé dans une première approche par des tractions surfaciques tangentielles réparties uniformément en fonction du rayon et un flux de chaleur uniforme sur la surface de contact.
## 3.1 Mise en données

La structure maillée (Figure 3(a)) possède un rayon de 35 mm et une épaisseur de 4 mm. Le problème compte 29105 degrés de liberté. Les conditions aux limites thermiques et mécaniques sont détaillées sur les figures 3(b) et 3(c) respectivement pour les faces supérieure et inférieure de la structure. Un flux de chaleur constant de 6 W/mm² et des tractions surfaciques tangentielles de 600 MPa sont imposés au niveau de la zone de frottement.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Maillage et conditions aux limites</figcaption></figure>

La température de transition fluide/solide est fixée à 100°C et la température extérieure est prescrite à 20°C. Les valeurs des paramètres matériau sont regroupées dans le tableau 1.

<table><caption>TABLE 1 – Paramètres matériau</caption><tr><td>Module d'Young <i>E</i></td><td>73000 MPa</td><td>Conductivité thermique <i>k</i></td><td>121 W/(m.K)</td></tr><tr><td>Coefficient de Poisson <i>ν</i></td><td>0.34</td><td>Masse volumique <i>ρ</i></td><td>2770 kg/m<sup>3</sup></td></tr><tr><td>Limite d'élasticité <i>σ</i><sub>y</sub></td><td>380 MPa</td><td>Capacité calorifique <i>C</i></td><td>875 J/(kg.K)</td></tr><tr><td>Module d'écrouissage <i>h</i></td><td>12000 MPa</td><td>Coefficient échange air <i>h</i><sub>air</sub></td><td>10 W/(m<sup>2</sup>.K)</td></tr><tr><td>Consistance (loi fluide) <i>K</i></td><td>10 MPa.s</td><td>Coefficient échange support <i>h</i><sub>supp</sub></td><td>350 W/(m<sup>2</sup>.K)</td></tr><tr><td>Sensibilité (loi fluide) <i>n</i></td><td>1</td><td></td><td></td></tr></table>

La limite d'élasticité et la consistance du fluide sont des fonctions décroissantes de la température.

## 3.2 Résultats

On présente sur la figure 4 l'évolution du champ de température sur la première partie du calcul. On observe bien l'apparition d'une zone chaude sous l'outil, qui commence d'abord à se propager sous forme d'une coupole hémisphérique puis de façon radiale sur la structure. En concordance avec le champ de température, la zone fluide apparaît peu après le début de la simulation au voisinage de la zone de contact outil/matière et progresse de façon hémisphérique au cours du temps. La figure 5 présente le champ de vitesse, on observe bien la transition solide/fluide entre les figures 5(a) et 5(b). Dès qu'une zone fluidifie, toute la déformation s'y concentre. Durant le procédé, l'outil créée des contraintes dans la zone où il appuie (Figure 6), mais ces contraintes chutent quand le solide devient fluide. Les contraintes sont donc maximales dans le solide au voisinage de l'interface fluide/solide.

## 4 Conclusions

Lors du procédé Friction Stir Spot Welding, deux phases sont générées, une pâteuse et l'autre demeurant solide. On a présenté un modèle basé sur un couplage fluide/solide développé avec l'élément P1+/P1, permettant de décrire correctement chacune des phases et, au travers de l'approche ALE, de suivre la matière dans la partie solide mais pas dans la partie fluide. Une simulation simplifiée de ce procédé a permis ici d'illustrer l'apparition puis la progression d'une phase fluide sous l'outil. Les vitesses atteintes dans la partie fluide lors de la simulation (Figure 5(c)) montrent que le modèle permet de simuler plusieurs tours de l'outil dans la matière.
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Évolution du champ de température.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Évolution de la norme de la vitesse matérielle.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Évolution de la contrainte équivalente de Mises.</figcaption></figure>

Références

[1] Arnold, D.N., Brezzi, F., Fortin, M. 1984 A stable finite element for the Stokes equations. *Calcolo*. **21** 337-344

[2] Feulvarch, E., Boitout, F., Bergheau, J.M. 2007 Friction Stir Welding : modélisation de l'écoulement de la matière pendant la phase de soudage. *European Journal of Computational Mechanics, Spécial*
« *Mécanique des fluides* » . **16** 865-887

[3] Heuzé, T., Leblond, J.B., Bergheau, J.M., Feulvarch, E. 2010 A finite element for laminar flow of incompressible fluids with inertia effects and thermomechanical coupling. *European Journal of Computational Mechanics*. **19** 293-304

[4] Sakano, R., Murakami, K., Yamashita, K., Hyoe, T., Fujimoto, M., Inuzuka, M., Nagao, U., Kashiki, H. 2001 Development of spot FSW robot system for automobile body members. *Third International Symposium of Friction Stir Welding*. Kobe, Japan.

[5] SYSWELD® 2007 User's Manual. ESI Group.

[6] Thomas, W.M., Nicholas, E.D., Smith, S.D. 1991 Proceedings of the TMS 2001 Aluminum Automotive and Joining Sessions, TMS. S.K. Das, J.G. Kaufman, T.J. Lienert (Eds.), *Aluminum 2001*.