# Meshless Simulation of Friction Stir Welding

I. Alfaro, L. Fratini, Elías Cueto, Francisco Chinesta, F. Micari

► To cite this version:

I. Alfaro, L. Fratini, Elías Cueto, Francisco Chinesta, F. Micari. Meshless Simulation of Friction Stir Welding. 9th International Conference on Numerical Methods in Industrial Froming Processes, 2007, Porto, Portugal. 10.1063/1.2740812. [hal-00109987](https://hal.science/hal-00109987v1)

HAL Id: hal-00109987

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

Submitted on 27 Aug 2020

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
# Meshless Simulation of Friction Stir Welding

I. Alfaro*, L. Fratini†, E. Cueto*, F. Chinesta** and F. Micari†

*Group of Structural Mechanics and Materials Modelling, Aragón Institute of Engineering Research (I3.4),
University of Zaragoza. María de Luna, 5. Campus Rio Ebro. E-50018 Zaragoza, Spain.

†Dept. of Technology, Production and Management Eng., University of Palermo, 90126 Palermo, Italy

**Laboratoire de Mécanique des Systèmes et des Procédés, LMSP UMR 8106 CNRS-ENSAM-ESEM. 105, Bvd. de l'Hôpital. F-75013 Paris, France.

**Abstract.** This paper encompasses our first efforts towards the numerical simulation of friction stir welding by employing a Lagrangian approach. To this end, we have employed a meshless method, namely the Natural Element Method (NEM). Friction Stir welding is a welding process where the union between the work pieces is achieved through the extremely high deformation imposed by a rotating pin, which moves between the two pieces. This extremely high strain is the main responsible of the difficulties associated with the numerical simulation of this forming process. Eulerian and Arbitrary Lagrangian-Eulerian (ALE) frameworks encounter difficulties in some aspects of the simulation. For instance, these approaches need additional techniques for the description of the boundary between materials, such as level sets, boundary markers or similar. In this paper we address the issue of employing a Lagrangian framework, which adequately describes the evolution in time of the interphase. The meshless character of the technique also ensures that no degeneracy on the accuracy is obtained due to mesh distortion. Some examples are presented that show the potential of the technique in simulating such a process.

**Keywords:** Natural Element Method, Friction Stir Welding, meshless.

## INTRODUCTION

Friction stir welding (FSW) is a process that, although in its development stage, has been successfully used to join pieces of materials with poor weldability. Briefly speaking, it achieves welding of the pieces by employing a rotating ping that provoques both extremely high plastic deformation and also a high heat generation. It is schematically represented in Fig. 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of the FSW process.</figcaption></figure>

In FSW, two pieces of sheet or thin plate are joined by inserting a specially designed rotating pin into the adjoining edges of the sheets to be welded and then moving it all along the seam. At first, the sheets or plates are abutted along edge to be welded and the rotating pin is sunken into the sheets/plates until the tool shoulder is in full contact with the sheets or plates surface. Once the pin is completely inserted, it is moved with a small nuting angle in the welding direction. Due to the advancing and rotating effect of the pin and shoulder of the tool

along the seam, an advancing side and a retreating side are formed and the softened and heated material flows around the pin to its backside where the material is consolidated to create a high-quality, solid-state weld.

In this papers we focus on the numerical simulation of such a process, rather than the experimental characterization. From the numerical point of view, such a process presents several challenging difficulties. First of all, the extremely large deformation present during the process, although very localized in a more or less region around the pin, always introduced numerical problems. Second, there is a strong coupling between these large deformation and heat generation, that in turn affects the behavior of the material.

Up to our knowledge, few numerical attempts have been made in order to simulate this process. In [4] a FE-based Lagrangian approach with intensive remeshing is employed. Similar approaches have been employed in [2] and [3]. Remeshing, however, is well-known as a potential source of numerical diffusion in the results.

In this paper we employ a somewhat different approach based on the use of meshless methods. Meshless methods allow for a Lagrangian description of the motion, while avoiding the need for remeshing. Thus, the nodes, that in our implementation transport all the variables linked to material's history, remain the same throughout the simulation. Mappings between old and new meshes are not necessary and hence the avoidance of numerical diffusion.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Delaunay triangulation and Voronoi diagram of a cloud of points.</figcaption></figure>

Among the many meshless methods available nowadays, we have chosen the Natural Element Method (NEM) [10] [6]. It possesses some noteworthy advantages over other meshless methods that will be described shortly.

## THE NATURAL ELEMENT METHOD

Consider a model composed by a cloud of points $N = \{n_1, n_2, \dots, n_m\} \subset \mathbb{R}^d$, for which there is a unique decomposition of the space into regions such that each point within these regions is closer to the node to which the region is associated than to any other in the cloud. This kind of space decomposition is called a Voronoi diagram (also Dirichlet tessellation) of the cloud of points and each Voronoi cell is formally defined as (see figure 2):

$$T_I = \{ \mathbf{x} \in \mathbb{R}^d : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) \quad \forall J \neq I \}, \quad (1)$$

where $d(\cdot, \cdot)$ is the Euclidean distance function.

The dual structure of the Voronoi diagram is the Delaunay triangulation, obtained by connecting nodes that share a common ($d - 1$)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there being some so-called *degenerate* cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in 2D, as depicted in Fig. 2 (right)). Another way to define the Delaunay triangulation of a set of nodes is by invoking the *empty circumcircle* property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called *natural neighbours* and hence the name of the technique.

Equivalently, the second-order Voronoi diagram of the cloud is defined as

$$T_{IJ} = \{ \mathbf{x} \in \mathbb{R}^d : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) < d(\mathbf{x}, \mathbf{x}_K) \quad \forall J \neq I \neq K \}. \quad (2)$$

The most extended natural neighbour interpolation method is the Sibson interpolant [8] [9]. Consider the introduction of the point $\mathbf{x}$ in the cloud of nodes. Due to this introduction, the Voronoi diagram will be altered, affecting the Voronoi cells of the natural neighbours of $\mathbf{x}$.

Sibson [8] defined the natural neighbour coordinates of a point $\mathbf{x}$ with respect to one of its neighbours $I$ as the ratio of the cell $T_I$ that is transferred to $T_X$ when adding $\mathbf{x}$ to the initial cloud of points to the total volume of $T_X$. In other words, if $\kappa(\mathbf{x})$ and $\kappa_I(\mathbf{x})$ are the Lebesgue measures of $T_X$ and $T_{XI}$ respectively, the natural neighbour coordinates of $\mathbf{x}$ with respect to the node $I$ is defined as

$$\phi_I(\mathbf{x}) = \frac{\kappa_I(\mathbf{x})}{\kappa(\mathbf{x})}. \quad (3)$$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Definition of the Natural Neighbour coordinates of a point $\mathbf{x}$.</figcaption></figure>

In Fig. 3 the shape function associated to node 1 at point $\mathbf{x}$ may be expressed as

$$\phi_1(\mathbf{x}) = \frac{A_{\bullet\bullet b\,f\,e}}{A_{\bullet\bullet b\,c\,d}}. \quad (4)$$

Sibson's interpolation scheme possesses the usual reproducing properties for this class of problems, i.e., verifies the *partition of unity* property (constant consistency), linear consistency (and therefore are suitable for the solution of second-order PDE). Other interesting properties such as the Kronecker delta property [10] and linear interpolation on the boundary [5] [11] are also verified by the NEM. This is especially important for problems involving friction or, in general, in which the compatibility along the boundary is important.

This kind of interpolation scheme is then used in a Galerkin framework, exactly like in the finite element method.

## CONSTITUTIVE MODELLING

FSW processes involve large deformation and high velocities of the rotating pin that make the elastic strains to be negligible. Although an elastic recovery exists, it is obvious negligible as a first approximation. The obvious advantage of this assumption is that the material can then be modelled as a non-newtonian (visco-plastic) fluid. This assumption is known as the *flow formulation*.
in the forming processes community [13]. Thus, the essential variables of the problem will be velocities and pressures, instead of displacements and pressures.

The deviatoric stresses will be, under this assumption,

$$ \mathbf{s} = 2\mu\mathbf{d}, \qquad (5) $$

being, as usual,

$$ \sigma = \mathbf{s} - p\mathbf{I} \qquad (6) $$

where $p = -tr(\sigma)/3$ and $\mathbf{I}$ stands for the second-order identity tensor. Obviously, in the most general case, the parameter $\mu$ will depend on both the level of strain (and hence the non-linear character of the behavior) and the temperature. To derive the expression of the parameter $\mu$ it is a common practice to write the strain rate tensor as emerging from a visco-plastic potential. Following Perzyna [7]

$$ \mathbf{d}^{vp} = \dot{\gamma} \frac{\partial Y(\sigma, \mathbf{q})}{\partial \sigma}, \qquad (7) $$

where $Y$ is the viscoplastic potential—usually coincident with the plastic criterion as has been considered here—, $\dot{\gamma}$ is a scalar function given by

$$ \dot{\gamma} = \frac{\langle g(Y(\sigma, \mathbf{q})) \rangle}{\eta} \quad \text{with} \quad \langle x \rangle = \frac{x + |x|}{2}, \qquad (8) $$

$\langle g \rangle$ is a monotonic function that takes zero value only if $Y(\sigma, \mathbf{q}) \le 0$, $\eta$ is a positive parameter often called viscosity and $\mathbf{q}$ represents the *hardening* parameters. In what follows we will avoid the use of the $vp$ superscript to indicate viscoplastic if there is no risk of confusion.

For metals, there exist well-defined plastic yield rules and for aluminium it is a common practice to employ a von Mises criterion:

$$ Y(\sigma, \mathbf{q}) = \bar{\sigma}(\mathbf{s}) - \sigma_y(\bar{\mathbf{d}}, T), \qquad (9) $$

where $\bar{\sigma} = \sqrt{\frac{3}{2}\mathbf{s} : \mathbf{s}} = \sqrt{3J_2}$ represents the *effective stress* and $\sigma_y$ represents the uniaxial yield stress. $\bar{\mathbf{d}}$ is the only internal variable in this model and is sometimes called *effective strain rate*:

$$ \bar{\mathbf{d}} = \sqrt{\frac{2}{3} \mathbf{d} : \mathbf{d}}. \qquad (10) $$

We have considered a yield stress given by the expression:

$$ \sigma_y = KT^A \bar{\mathbf{d}}^B \bar{\mathbf{\epsilon}}^C \qquad (11) $$

where $K = 2.69E10$, $A = -3.3155$, $B = 0.1324$ and $C$ is usually taken as 0.0192. In this work we have neglected the influence of strain, thus taking $C = 0$ instead. This law has been employed in previous works on this topic, see [2].

If we combine now the general form of the strain rate tensor given in Eq. (7), with Eq. (9), we arrive to

$$ \mathbf{d} = \dot{\gamma} \frac{3\mathbf{s}}{2\bar{\sigma}}. \qquad (12) $$

It is immediate now, by combining Eq. (10) and the definition of effective stress, $\bar{\sigma}$, to prove that $\dot{\gamma}$ is precisely the effective strain rate:

$$ \bar{\mathbf{d}} = \frac{\dot{\gamma}}{\bar{\sigma}} \sqrt{\frac{3}{2} \mathbf{s} : \mathbf{s}} = \dot{\gamma}. \qquad (13) $$

On the other hand, and by following the Perzyna-like model employed in Eqs. (7) and (8) and taking $g(f) = f$, we arrive to a relationship between equivalent stress and equivalent strain rate:

$$ \bar{\mathbf{d}} = \dot{\gamma} = \frac{\langle \bar{\sigma} - \sigma_y \rangle}{\eta} \Rightarrow \bar{\sigma} = \eta \bar{\mathbf{d}} + \sigma_y \text{ if } \bar{\sigma} \ge \sigma_y \qquad (14) $$

that, introduced in Eq. (12), accounting Eq. (13), gives the following visco-plastic constitutive equation:

$$ \mathbf{s} = 2 \frac{\eta \bar{\mathbf{d}} + \sigma_y(\bar{\mathbf{d}})}{3\bar{\mathbf{d}}} \mathbf{d}. \qquad (15) $$

Note that, depending on the $\eta$ value, the return to the yield surface is done with different velocity. Since it is common to describe aluminium behaviour as rigid-plastic (rather than viscoplastic) we employ null viscosity $\eta$, so as to enforce $Y = \bar{\sigma} - \sigma_y = 0$, leading to

$$ \mathbf{s} = \frac{2\sigma_y}{3\bar{\mathbf{d}}} \mathbf{d}. \qquad (16) $$

Finally, the constitutive equation, accounting the incompressibility of plastic flow results:

$$ \sigma = 2\mu\mathbf{d} - p\mathbf{I}, \text{ with } \mu = \frac{\sigma_y}{3\bar{\mathbf{d}}}. \qquad (17) $$

## Governing equations

We consider the balance of momentum equations, without inertia and mass terms for the plates being welded—it is obviously not the case for the pin, which is considered as perfectly rigid as a first approximation—

$$ \nabla \cdot \sigma = 0, \qquad (18) $$

and the assumed incompressibility of a von Mises-like flow:

$$ \nabla \cdot \mathbf{v} = 0, \qquad (19) $$

where $\mathbf{v}$ represents the velocity field. The stress-strain rate relationship is given by Eq. (16). Velocities are interpolated by means of the functions 3, while pressures
are considered constant over the whole Voronoi cell associated to each node.

This rigid-plastic material is coupled with the following heat transfer equations:

$$ \nabla(k\nabla T) + \dot{r} - (\rho c_p \dot{T}) = 0 \quad (20) $$

where k denotes the thermal conductivity, $\dot{r}$ the heat generation rate, $\rho$ the specific density and $c_p$ the specific heat of the metal. The rate of heat generation in the aluminium billet due to plastic deformation is given by

$$ \dot{r} = \beta \sigma : \mathbf{d} \quad (21) $$

where $\beta$ represents the fraction of mechanical energy transformed to heat and is assumed to be 0.9 [12].

The coupling has been made through a block-iterative semi-implicit method, together with a fixed-point algorithm to treat the non-linear coupling. This strategy has also been employed successfully in previous works of the authors [1].

## NUMERICAL RESULTS

We consider here a very simplified two-dimensional, plane stress, model of the union of two plates. The model is composed by 240 nodes whose Delaunay triangulation is re-computed at each time-step. This is very fast process, which actually consumes only a few time, if compared with the equilibrium iterations. The pin is assumed to rotate at 1000 rev/min and to advance at 1mm/s. Despite that the usual friction coefficient between the pin and the sheet has been established to be around 0.5 (0.46 in [2]) we have considered here perfect adhesion between them. This will overestimate the strain produced by the pin rotation. The time increment chosen was $10^{-4}$ seconds.

In Figs. 4-7 the evolution of the temperature is shown at different time steps. In Figs. 8-11, respectively, the equivalent strain rate is depicted. In all cases results show an overall good qualitative agreement with those in [2], performed by Finite Element simulation.

However, the most important advantage of the presented technique is the possibility of tracking the materials particles throughout the process. We have considered two plates of different materials. In Fig. 12 the initial configuration of the nodes is depicted. In it, it can be shown that nodes of each plates have been labelled, by painting them in blue and red, respectively.

The simulation has been run for approximately 4000 time steps. Configuration after 4000 time steps is shown in Fig. 13.

This tracking can be performed without remeshing despite the high level of distortion of the resulting triangulation because triangles do not constitute the support of

<figure><img src="image_5.png" /><figcaption>Figure 4. Evolution of the temperature field during the beginning of the rotation. 10th time step.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of the temperature field during the beginning of the rotation. 50th time step.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Evolution of the temperature field during the beginning of the rotation. 70th time step.</figcaption></figure>

NEM shape functions. In Fig. 14 the resulting triangulation at the 4000th time step is shown. In this way, the numerical diffusion provoked by continuous remeshing is avoided. In addition, if some kind of nodal integration is performed, all history variables can be linked to the
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Evolution of the temperature field during the beginning of the rotation. 90th time step.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 10.** Evolution of the equivalent strain rate during the beginning of the rotation. 70th time step.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Evolution of the equivalent strain rate during the beginning of the rotation. 10th time step.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 11.** Evolution of the equivalent strain rate during the beginning of the rotation. 90th time step.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Evolution of the equivalent strain rate during the beginning of the rotation. 50th time step.</figcaption></figure>

nodes, thus avoiding any projection procedures.

## CONCLUSIONS

In this paper we have presented the results of the first attempt —up to our knowledge— of applying meshless methods to the simulation of Friction Stir Welding

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Initial configuration of the simulation with two materials.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Evolution of the phases through the welding process. 4000th time step</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Delaunay triangulation at the 4000th time step of the simulation.</figcaption></figure>

processes. The use of meshless methods (in particular, the Natural Element Method) was motivated trying to avoid extensive remeshing associated with the large deformations present in this kind of processes.

Although the presented results can be considered as very preliminary, and further refinement of the models (especially relative to contact and friction) is needed, we can conclude that the Natural Element Method constitutes a valuable tool for the simulation of such a complex forming process. In particular, it has been shown how a large number of time steps (up to four thousands) have been accomplished maintaining the initial set of nodes. Despite the high distortion of the triangulation, good qualitative agreement with previous results has been found.

Note, however, that the true interest of the proposed technique relies in its extension to three-dimensional settings, where the true complexity of the process should be analyzed. This constitutes the ongoing work of the authors.

## REFERENCES

1. I. Alfaro, D. Bel, E. Cueto, M. Doblaré, and F. Chinesta. Three-dimensional simulation of aluminium extrusion by the alpha-shape based natural element method. *Computer Methods in Applied Mechanics and Engineering*, 195(33-36):4269-4286, **2006**.

2. G. Buffa, J. Hu, R. Shivpuri, and L. Fratini. A continuum based fem model for friction stir welding-model development. *Materials Science and Engineering A*, 419:389-396, **2006**.

3. G. Buffa, J. Hu, R. Shivpuri, and L. Fratini. Design of the friction stir welding tool using the continuum based FEM model. *Materials Science and Engineering A*, 419:381-388, **2006**.

4. J.-L. Chenot and E. Massoni. Finite element modelling and control of new metal forming processes. *International Journal of Machine Tools and Manufacture*, 46:1194-1200, **2006**.

5. E. Cueto, M. Doblaré, and L. Gracia. Imposing essential boundary conditions in the Natural Element Method by means of density-scaled $\alpha$-shapes. *International Journal for Numerical Methods in Engineering*, 49-4:519-546, **2000**.

6. E. Cueto, N. Sukumar, B. Calvo, M. A. Martínez, J. Cegoñino, and M. Doblaré. Overview and recent advances in Natural Neighbour Galerkin methods. *Archives of Computational Methods in Engineering*, 10(4):307-384, **2003**.

7. P. Perzyna. Fundamental problems in visco-plasticity. In *Recent advances in Applied Mechanics*. Academic Press, New York, 1966.

8. R. Sibson. A Vector Identity for the Dirichlet Tesselation. *Mathematical Proceedings of the Cambridge Philosophical Society*, 87:151-155, 1980.

9. R. Sibson. A brief description of natural neighbour interpolation. In *Interpreting Multivariate Data. V. Barnett (Editor)*, pages 21-36. John Wiley, 1981.

10. N. Sukumar, B. Moran, and T. Belytschko. The Natural Element Method in Solid Mechanics. *International Journal for Numerical Methods in Engineering*, 43(5):839-887, 1998.

11. J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta. A new extension of the Natural Element method for non-convex and discontinuous problems: the Constrained Natural Element method. *International Journal for Numerical Methods in Engineering*, 60(8):1452-1474, **2004**.

12. J. Zhou, L. Li, and J. Duszczyk. 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach. *Journal of Materials Processing Technology*, 134:383-397, **2003**.

13. O. C. Zienkiewicz and P. N. Godbolet. Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes. *International Journal for Numerical Methods in Engineering*, 8:3-16, 1974.