# Formulation instationnaire d'un élément P1+/P1 pour la modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température

Thomas Heuzé, Jean-Baptiste Leblond, Jean-Michel Bergheau, Éric Feulvarch

► To cite this version:

Thomas Heuzé, Jean-Baptiste Leblond, Jean-Michel Bergheau, Éric Feulvarch. Formulation instationnaire d'un élément P1+/P1 pour la modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413223

HAL Id: hal-01413223

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

Submitted on 9 Dec 2016

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
# Formulation instationnaire d'un élément P1+/P1 pour la modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température

T. Heuzé<fnref n="1" />, J.B. Leblond<fnref n="1" />, J.M. Bergheau<fnref n="2" />, É. Feulvarch<fnref n="2" />

**Résumé** — Le procédé de soudage ponctuel par friction et malaxage fait intervenir de grandes déformations au voisinage de l'outil. La simulation de ce procédé doit rendre compte d'une phase pâteuse dans laquelle la matière est malaxée, et d'une phase demeurant solide. Une approche de type eulérienne lagrangienne arbitraire combinée à des comportements respectivement fluide et solide dans chacune de ces phases peut permettre de simuler de nombreux tours de l'outil dans la matière tout en suivant les bords des tôles. Le travail présenté ici se focalise sur une première étape de cette étude, le développement d'une formulation mixte température/vitesse/pression d'un élément fini fluide P1+/P1 dans le cas instationnaire.

**Mots clés** — Élément fini P1+/P1, Friction Stir Spot Welding, Instationnarité.

## 1 Introduction

Le procédé de soudage ponctuel par friction et malaxage (Friction Stir Spot Welding [8]) est un procédé dérivé du soudage par friction et malaxage (Friction Stir Welding [10]) consistant à créer un joint de soudure « ponctuel » entre deux tôles 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. Le mélange créé assure alors la solidarité des parties à joindre.

Les approches purement lagrangienne et eulérienne ne sont pas appropriées pour modéliser numériquement ce procédé, l'une nécessitant de fréquents remaillages dans le cas de grandes déformations et l'autre ne permettant pas de suivre les bords des tôles pendant le procédé. 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, peut permettre de simuler de nombreux tours de l'outil dans la matière tout en suivant les bords des tôles. En plus du traitement numérique différent des deux zones, on en profite pour distinguer les modèles de comportement dans chacune des zones en y intégrant respectivement des descriptions de type fluide visqueux incompressible et solide.

Dans une première étape de cette étude, on met en place l'élément de la partie fluide. L'élément fini P1+/P1 [1], basé sur une formulation mixte, est utilisé pour cette simulation. Une version déjà

<footnote n="1">¹ Institut Jean Le Rond D'Alembert, CNRS UMR 7190, Université Pierre et Marie Curie
4, Place Jussieu, tour 55-65, 75252 Paris Cedex 05
thomas.heuze@etu.upmc.fr, leblond@lmm.jussieu.fr</footnote>

<footnote n="2">² LTDS, UMR 5513 CNRS/ECL/ENISE
58, rue Jean Parot 42023 Saint Étienne Cedex 2
bergheau@enise.fr, feulvarch@enise.fr</footnote>
implémentée dans le code SYSWELD [9, 5] avec une formulation température/vitesse/pression pour un état stationnaire, est ici étendue au cas instationnaire. Ce travail présente la formulation de cet élément et sa validation sur un cas test, le viscosimètre de Couette.

## 2 Formulation mixte instationnaire température/vitesse/pression pour l'élément fluide P1+/P1

### 2.1 L'élément P1+/P1

La modélisation de la zone pâteuse comme un fluide visqueux incompressible fait intervenir la pression hydrostatique comme inconnue à part entière. Du point de vue de la simulation numérique, une discrétisation particulière doit être mise en oeuvre pour éviter des phénomènes numériques indésirables [7]. Le tétraèdre P1+/P1 [1] est un élément fini linéaire (Figure 1) assurant, dans la version proposée, la continuité des champs de température, vitesse et pression. L'exposant « + » réfère au noeud bulle ajouté au centre de l'élément, celui-ci permet d'enrichir le champ de vitesse et donc de satisfaire correctement les équations du mouvement et l'incompressibilité.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Tétraèdre P1+/P1 de référence.</figcaption></figure>

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)\lambda \quad 1 \le p \le 4 \quad (1) $$

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

### 2.2 Forme faible

Pour le procédé Friction Stir Spot Welding [8], un fort couplage thermomécanique doit être pris en compte. Le problème à résoudre fait alors intervenir les équations d'énergie, de mouvement et la liaison interne. Dans le cas d'une formulation mixte, le problème s'écrit sous forme faible :

$$ (W) \left\{ \begin{array}{l} \text{Trouver}(T, \mathbf{v}, p) \in (\mathcal{T}_{ad} \times \mathcal{V}_{ad} \times \mathcal{P}_{ad}), \forall t \in ]0, T], \text{tels que } \forall (\theta, \mathbf{w}, q) \in (\mathcal{T}_{ad}^0 \times \mathcal{V}_{ad}^0 \times \mathcal{P}_{ad}), \\ \left| \begin{array}{l} (\rho C \dot{T}, \theta) + A(T; T, \theta) + B(\mathbf{v}; T, \theta) = L(\theta) + (\mathcal{D}(\mathbf{v}), \theta) \\ (\rho \dot{\mathbf{v}}, \mathbf{w}) + a(\mathbf{v}, T; \mathbf{v}, \mathbf{w}) + b(\mathbf{v}; \mathbf{v}, \mathbf{w}) - (p, \text{div} \mathbf{w}) = l(\mathbf{w}) \\ -(q, \text{div} \mathbf{v}) = 0 \end{array} \right. \\ T(\mathbf{x}, t = 0) = T_0(\mathbf{x}) \\ \mathbf{v}(\mathbf{x}, t = 0) = \mathbf{v}_0(\mathbf{x}) \end{array} \right. \quad (2) $$
$$
\begin{align}
\text{avec } \bar{a}(\mathbf{v}, T; \mathbf{v}, \mathbf{w}) &= \int_{\Omega} 2\mu(T, \mathbf{v}) D_{ij}(\mathbf{v}) D_{ij}(\mathbf{w}) d\Omega; \quad A(T; T, \theta) = \int_{\Omega} k_{ij}(T) T_{,j} \theta_{,i} d\Omega \nonumber \\
b(\mathbf{v}; \mathbf{v}, \mathbf{w}) &= \int_{\Omega} \rho v_j v_{i,j} w_i d\Omega; \quad B(\mathbf{v}; T, \theta) = \int_{\Omega} \rho C v_i T_{,i} \theta d\Omega \tag{3} \\
l(\mathbf{w}) &= \int_{\Omega} f_i w_i d\Omega + \int_{\partial\Omega_T} T_i^d w_i dS; \quad L(\theta) = \int_{\Omega} r \theta d\Omega + \int_{\partial\Omega_q} \phi^d \theta dS \nonumber
\end{align}
$$

La notation (., .) désigne le produit des quantités dans L², D est la dissipation générée et φ^d est le flux thermique imposé. On note qu'une régularité plus forte est requise pour les champs de température et de vitesse que pour le champ de pression.

2.3 Formulation éléments finis

Le problème (W) est discrétisé par éléments finis. Un algorithme aux différences finies de type
Euler implicite est utilisé pour la résolution, l'intervalle de temps considéré est divisé en N pas de
temps Δt_n, n = 1 à N, et T = ΣN/2 Δt_n. L'algorithme de type implicite utilisé ici et la non linéarité
de ce problème nécessitent une résolution itérative. 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}_{n+1}^{(k)} \dot{\mathbf{q}}_{n+1}^{(k+1)} + \mathbf{K}_{n+1}^{(k)} \delta \mathbf{q}^{(k)} = \mathbf{R}_{n+1}^{(k)} \quad (4)
$$

où **M**, **K** et **R** sont respectivement les matrices de masse généralisée, de raideur tangente géné-
ralisée et le vecteur résidu, évalués à l'itération *k* du calcul. Le vecteur des degrés de liberté du
système **q** est défini par :

$$
\mathbf{q}^T = \{ \mathbf{T} \ \mathbf{v} \ \mathbf{p} \ \lambda \} \tag{5}
$$

et son incrément $\delta \mathbf{q}^{(k)}$ par $\delta \mathbf{q}^{(k)} = \mathbf{q}_{n+1}^{(k+1)} - \mathbf{q}_{n+1}^{(k)}$. Le système d'équations (4) est schématiquement composé de la manière suivante :

$$
\begin{bmatrix} \mathbf{C} & 0 & 0 & 0 \\ 0 & \mathbf{M}_{vv} & 0 & \mathbf{M}_{vb} \\ 0 & 0 & 0 & 0 \\ 0 & \mathbf{M}_{vb}^T & 0 & \mathbf{M}_{bb} \end{bmatrix}_{n+1}^{(k)} \begin{Bmatrix} \dot{\mathbf{T}} \\ \dot{\mathbf{v}} \\ \dot{\mathbf{p}} \\ \dot{\lambda} \end{Bmatrix}_{n+1}^{(k+1)} + \begin{bmatrix} \mathbf{K}_{TT} & \mathbf{K}_{T\nu} & 0 & \mathbf{K}_{Tb} \\ \mathbf{K}_{\nu T} & \mathbf{K}_{\nu\nu} & \mathbf{K}_{\nu p} & \mathbf{K}_{\nu b} \\ 0 & \mathbf{K}_{p\nu} & 0 & \mathbf{K}_{pb} \\ \mathbf{K}_{bT} & \mathbf{K}_{bv} & \mathbf{K}_{bp} & \mathbf{K}_{bb} \end{bmatrix}_{n+1}^{(k)} \begin{Bmatrix} \delta\mathbf{T} \\ \delta\mathbf{v} \\ \delta\mathbf{p} \\ \delta\lambda \end{Bmatrix}^{(k)} = \begin{Bmatrix} \mathbf{R}_T \\ \mathbf{R}_\nu \\ \mathbf{R}_p \\ \mathbf{R}_\lambda \end{Bmatrix}_{n+1}^{(k)} \quad (6)
$$

où **K** est assemblée à partir des contributions des sous-matrices **K**<sub>xy</sub> = ∂**R**<sub>x</sub>/∂**y**, (x, y) = (T, v, p, λ). Les termes convectifs mécaniques et thermiques sont contenus dans le second membre, ils nécessitent par ailleurs un traitement particulier [5, 2, 3] pour assurer la stabilité des solutions numériques. De plus, l'hypothèse est faite que l'accélération exprimée au noeud bulle est environ égale à la moyenne de celle des noeuds extrêmes, il en résulte que λ̃ ≃ 0. Enfin, pour un gain de temps lors de la résolution, le noeud bulle est condensé dans les équations relatives aux noeuds extrêmes. La résolution de ce système s'effectue alors en résolvant un système linéaire du type :

$$
(\mathbf{K}^*)_{n+1}^{(k)} \mathbf{q}_{n+1}^{(k+1)} = (\mathbf{R}^*)_{n+1}^{(k)} \quad (7)
$$

3 Validation de l'élément P1+/P1 en instationnaire : application au viscosimètre de Couette

La validation de la formulation de cet élément est effectuée sur le problème du viscosimètre de
Couette. Le problème défini est présenté sur la figure 2. Deux solutions analytiques instationnaires
sont présentées respectivement en mécanique instationnaire et pour un couplage (thermique installationnaire - mécanique stationnaire). Les solutions sont obtenues en prenant des conditions aux limites et initiales favorisant une expression des champs solution la plus simple possible.

## 3.1 Résolution mécanique instationnaire

Les hypothèses de la résolution mécanique instationnaire sont les suivantes :

- le problème est plan, les champs sont invariants par rapport à la coordonnée z,

- le fluide est non pesant,

- le problème présente une symétrie de révolution,

- le fluide est considéré comme Newtonien et l'écoulement comme incompressible,

- le cylindre extérieur est fixé, le mouvement du cylindre intérieur est imposé.

De l'équation tangentielle du mouvement découle le calcul du champ de vitesse. Par une méthode de séparation des variables, on obtient une équation de Bessel [4] dont la solution est une combinaison linéaire de fonctions de Bessel de première et seconde espèce du premier ordre dans notre cas. Pilotant la vitesse du cylindre intérieur et en choisissant une condition initiale non nulle convenable, il vient le champ de vitesse transitoire :

$$v(r,t) = v_{\theta}(r,t)\mathbf{e}_{\theta} = \Omega_{10}a \left( \frac{J_1(r)Y_1(b) - J_1(b)Y_1(r)}{J_1(a)Y_1(b) - J_1(b)Y_1(a)} \right) \exp\left(-\frac{\mu}{\rho}t\right)\mathbf{e}_{\theta} \quad (8)$$

où $a$ et $b$ sont les rayons des cylindres intérieur et extérieur, $\rho$ et $\mu$ la masse volumique et la viscosité dynamique du fluide, $\Omega_{10}$ la vitesse de rotation initiale du cylindre intérieur, $J_1$ et $Y_1$ les fonctions de Bessel de première et seconde espèce d'ordre 1.

## 3.2 Résolution thermique instationnaire / mécanique stationnaire

L'écriture du bilan d'énergie interne fait ici apparaître la dissipation comme terme source. Ce terme est calculé à partir d'un problème choisi de manière à ce que la résolution mécanique soit indépendante de la résolution thermique et puisse donc être effectué à priori. Au moyen d'un changement de variable approprié, il est possible de se ramener à une équation sans terme source sur une nouvelle variable $\bar{T}$, la résolution s'effectue alors de la même façon que précédemment. Des conditions homogènes de Dirichlet et Neumann sont respectivement imposées en $r = a$ et $r = b$ sur la variable $\bar{T}$. La solution thermique instationnaire s'exprime alors ainsi :

$$ 
\begin{gathered}
T(r,t) = \left( J_0(\lambda_1 r) + \frac{Y_0(\lambda_1 r)}{\eta} \right) \exp(-\alpha t) - \frac{A}{4kr^2} \\
\text{avec } \eta = -\frac{Y_0(\lambda_1 a)}{J_0(\lambda_1 a)}, \quad \alpha = \kappa\lambda_1^2, \quad \kappa = \frac{k}{\rho C}, \quad A = 4\mu\Omega_1^2 \frac{(ab)^2}{(a^2 - b^2)^2} \\
\text{et } \left| \begin{array}{l} T(a,t) = -\frac{A}{4ka^2} \\ \frac{\partial T}{\partial r}(b,t) = \frac{A}{2kb^3} \end{array} \right.
\end{gathered}
\quad (9) $$

où $k$, $\rho$, $C$ sont la conductivité thermique, la masse volumique et la capacité thermique du fluide, $\Omega_1$ est la vitesse de rotation du cylindre intérieur et $\lambda_1$ est la première racine de l'équation définie par les conditions aux limites. On observe que les conditions homogènes appliquées sur $\bar{T}$ conduisent à des conditions non homogènes par changement inverse, la température imposée en $r = a$ est alors négative.
## 3.3 Comparaison des résultats numériques/analytiques

Le domaine fluide est discrétisé avec l'élément P1+/P1, un maillage est présenté sur la figure 3. Les rayons intérieur et extérieur sont fixés tels que $a = 0.1\text{m}$ et $b = 1\text{m}$. Par convenance, les paramètres matériau du fluide sont tous prescrits à l'unité.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Viscosimètre de Couette</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3 – Maillage du viscosimètre de Couette</figcaption></figure>

### 3.3.1 Mécanique instationnaire

Les valeurs des conditions initiales analytiques sont spécifiées en chaque noeud dans la mise en données, et la vitesse des noeuds du cylindre intérieur est imposée suivant une exponentielle décroissante. La thermique a elle été bloquée, c'est-à-dire qu'une température nulle est imposée à tous les noeuds. La figure 4 présente la comparaison analytique/numérique pour la solution mécanique transitoire où les valeurs à quelques pas de temps sont présentées. Le champ de vitesse est tracé sur un rayon du viscosimètre, les points repérés sur la courbe correspondent aux solutions éléments finis extraites aux noeuds du maillage sur une ligne radiale visible sur la figure 3. On observe une très bonne corrélation sur les instants intermédiaires entre les deux méthodes.

### 3.3.2 Thermique instationnaire - mécanique stationnaire

Des conditions initiales thermiques et mécaniques issues des solutions analytiques sont spécifiées dans la mise en donnée. Par la suite, le cylindre intérieur est piloté à vitesse de rotation constante $\Omega_1 = 10 \text{rd/s}$, et les conditions de température et flux (9) sont prescrites sur les noeuds des cylindres intérieur et extérieur. La figure 5 présente la comparaison analytique/numérique pour la solution thermique transitoire. Le champ de température présenté est tracé sur une ligne radiale du viscosimètre. On observe aussi pour la température une bonne corrélation entre les deux solutions au cours des pas de temps du calcul.

## 4 Conclusion

Lors du process Friction Stir Spot Welding, deux phases sont générées, une pâteuse et l'autre demeurant solide. Dans ce travail, un premier outil pour la simulation de la phase pâteuse a été mis en place. L'élément P1+/P1 développé dans le cadre d'une formulation mixte température/vitesse/pression est ici étendu au cas instationnaire. La validation de cet élément est effectuée sur le problème du viscosimètre de Couette, deux solutions analytiques transitoires sont construites
<figure><img src="image_5.png" /><figcaption>Figure 4 – Champ de vitesse transitoire.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5 – Champ de température transitoire.</figcaption></figure>

pour tester les comportements mécanique et thermique de l'élément. La comparaison des solutions analytiques et numériques fait apparaître un comportement sain de l'élément sur ces problèmes transitoires.

Références

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

[2] T. Belytschko, W.-K. Liu, B. Moran. *Nonlinear Finite Elements for Continua and Structures*, Wiley, 2000.

[3] Jean-Michel Bergheau, Roland Fortunier. *Simulation numérique des transferts thermiques par éléments finis*, Hermès, 2004.

[4] H.-S. Carslaw, J.-C. Jaeger. *Conduction of heat in solids*, Oxford, Clarendon press – cop, Second edition, 1959.

[5] Éric Feulvarch. *Modélisation numérique du soudage par friction-malaxage (Friction Stir Welding)*, thèse de l'Université Jean Monnet de Saint Étienne, 2005.

[6] Simon Guerdoux. *Numerical simulation of the Friction Stir Process*, thèse de l'École des Mines de Paris, 2007.

[7] Thomas J.R. Hughes. *The Finite Element Method : Linear Static and Dynamic Finite Element Analysis*, Dover, 1987.

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

[9] SYSWELD®. *User's Manual*, ESI Group, 2007.

[10] W.-M. Thomas, E.-D. Nicholas, J.-C. Nedd Ham, M.-G. Murch, P. Templesmith, C.-J. Dawes. Friction stir butt welding. *International patent application no. PCT/GB92/02203 and GBpatent application no. 9125978.8*, 6 December 1991.