Meshless Simulation of Friction Stir Welding

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). Fric­ tion Stir welding is a welding process where the union between the work pieces is achieved through the extremely high defor­ mation 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.


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.Briefl y speaking, it achieves welding of the pieces by employ ing a rotating ping that provoques both extremely high plastic deformation and also a high heat generation.It is schematically represented in Fig. 1.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.

CD Piece 1
In this papers we focus on the numerical simulation of such a process, rather than the experimental characteri zation.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 deforma tion and heat generation, that in tum aff ects 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 ap proach based on the use of meshless methods.Mesh less 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.Among the many meshless methods available nowa days, we have chosen the Natural Element Method (NEM) [10] [6].It possesses some noteworthy advan tages over other meshless methods that will be described shortly.

THE NATURAL ELEMENT METHOD
Consider a model composed by a cloud of points N = { n1,n2, ... ,nm} c JR. d , for which there is a unique decom position 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): 0={xElR d :d(x,x1)<d(x,xJ)VJ#I}, (1)   where d( •, •) is the Euclidean distance function.
The dual structure of the Voronoi diagram is the De launay triangulation, obtained by connecting nodes that share a common (d -I)-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 triangu lations (consider, for example, the case of triangulating a square in 2D, as depicted in Fig. 2  Equivalently, the second-order Voronoi diagram of the cloud is defined as The most extended natural neighbour inteipolation method is the Sibson interpolant [8] [9].Consider the introduction of the point x in the cloud of nodes.Due to this introduction, the Voronoi diagram will be altered, aff ecting the Voronoi cells of the natural neighbours of x. Sibson [8] defined the natural neighbour coordinates of a point x with respect to one of its neighbours I as the ratio of the cell TJ that is transferred to Tx when adding x to the initial cloud of points to the total volume of Tx.In other words, if K ( x) and K1 ( x) are the Lebesgue measures of Tx and TxI respectively, the natural neighbour coordinates of x with respect to the node I is defined as (3) In Fig. 3 the shape function associated to node 1 at point x may be expressed as Sibson's interpolation scheme possesses the usual re producing properties for this class of problems, i.e., veri fies the partition of unity property (constant consistency), linear consistency (and therefore are suitable for the so lution of second-order PDE).Other interesting properties such as the Kronecker delta property [10] and linear in terpolation on the boundary [5] [11] are also verified by the NEM.This is especially important for problems in volving 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 iINolve large deformation and high ve locities of the rotating pin that make the elastic strains to be negligible.Although an elastic recovery exists, it is obvious negligible as a fi rst approximation.The ob vious 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 fo rmulation in the forming processes community [13].Thus, the es sential variables of the problem will be velocities and pressures, instead of displacements and pressures.
The deviatoric stresses will be, under this assumption, s = 2µd, being, as usual, where p = -tr( <J ) /3 and I stands for the second-order identity tensor.Obviously, in the most general case, the parameterµ 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 µ it is a common practice to write the strain rate tensor as emerging from a visco-plastic potential.Following Perzyna [7] dvp -. dY ( <J , q ) -r d<J , where Y is the viscoplastic potential -usually coincident with the plastic criterion as has been considered here-, r is a scalar function given by i' = (g( Y ( ;, q ))) with ( x (g) is a monotonic function that takes zero value only if Y ( CJ, q ) ::; 0, 17 is a positive parameter often called viscosity and 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 ( CJ, q ) = O'( s ) -<J y (d, T), (9) where O' = � = VJJ2 represents the effective stress and <J y represents the uniaxial yield stress.d is the only internal variable in this model and is sometimes called effective strain rate: (10) We have considered a yield stress given by the expression: (11) where K = 2. 69£10, 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 It is immediate now, by combining Eq. ( 10) and the def inition of effective stress, O', to prove that r is precisely the effective strain rate: -. � d= �y2s:s=Y. (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: that, introduced in Eq. ( 12), accounting Eq. ( 13), gives the following visco-plastic constitutive equation: (15) Note that, depending on the 17 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 viscos  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-Y'.(J = 0, (18) and the assumed incompressibility of a von Mises-like fl ow: where v represents the velocity field.The stress-strain rate relationship is given by Eq. ( 16).Velocities are in terpolated by means of the functions 3, while pressures are considered constant over the whole Voronoi cell as sociated to each node.This rigid-plastic material is coupled with the follow ing heat transfer equations: where k denotes the thermal conductivity, r the heat generation rate, p 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 r=f3cr:d where f3 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 algo rithm to treat the non-linear coupling.This strategy has also been employed successfully in previous works of the authors [l].

NUMERICA L 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 triangula tion 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 Imm/ s.Despite that the usual friction coefficient be tween 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 pre sented technique is the possibility of tracking the ma terials particles throughout the process.We have consid ered two plates of different materials.In Fig. 12 the ini tial 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 inFig.13.This tracking can be performed without remeshing de spite the high level of distortion of the resulting triangu lation because triangles do not constitute the support of  NEM shape functions.In Fig. 14 the resulting triangu lation 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   processes.The use of meshless methods (in particular, the Natural Element Method) was motivated trying to avoid extensive remeshing associated with the large de formations 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 consti tutes a valuable tool for the simulation of such a com plex forming process.In particular, it has been shown how a large number of time steps (up to four thou sands) have been accomplished maintaining the initial set of nodes.Despite the high distortion of the triangu lation, 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.

Figure 1 .
Figure 1.Schematic representation of the FSW process.

Figure 2 .
Figure 2. Delaunay triangulation and Voronoi diagram of a cloud of points.
(right)).Another way to define the Delaunay triangulation of a set of nodes is by iINoking 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.

Figure 3 .
Figure 3. Definition of the Natural Neighbour coordinates of a point x.

Figure 4 .
Figure 4. Evolution of the temperature field during the be ginning of the rotation.lOth time step.

Figure 5 .Figure 6 .
Figure 5. Evolution of the temperature field during the be ginning of the rotation.50th time step. x

x 50 Figure 13 .Figure 14 .
Figure 13.Evolution of the phases through the welding process.4000th time step