An Implicit Algorithm Based on Continuous Moving Least Square to Simulate Material Mixing in Friction Stir Welding Process

An implicit iterative algorithm, based on the continuous moving least square (CMLS), is developed to simulate material mixing in Friction Stir Welding (FSW) process. Strong formulation is chosen for the modeling of the mechanical problem in Lagrangian framework to avoid the drawback of numerical integration.This algorithm is well adapted to large deformations in the mixing zone in the neighborhood of the welding tool. We limit ourselves to bidimensional viscoplastic problem to show the performance of the proposed implicit algorithm. The results show that the proposed algorithm can be employed to simulate FSW.


Introduction
Friction stir welding (FSW) [1] is a solid state welding process used to join two pieces (plates, tubes, spheres, etc.) of aluminum alloys.Joining two workpieces by FSW consists in heat generation mainly due to the shoulder and material mixing thanks to the pin that provokes both extremely high plastic deformation and also a high heat generation.
Numerical simulation of friction stir welding process is important in industry to control the product quality and optimize the fabrication cost and avoid the traditional tests.It has been investigated by several authors considering thermal or thermomechanical framework [2][3][4][5][6][7][8].The difficulty in modeling this welding type resides to large deformations generated by the material mixing.Different formulations have been proposed in Eulerian, Lagrangian, or Arbitrary Lagrangian Eulerian frameworks.Eulerian formulation is appropriate to describe material flow for large deformations thanks to fixed mesh grid but this formulation cannot be used for problems involving free surfaces.Lagrangian formulation is well adapted for history dependence in solids and for simulation of material flow with free surfaces but mesh distortions require special procedure as remeshing.The ALE formulation has been proposed to benefit from both advantages of Eulerian and Lagrangian formulations [3].In this technique, one has to describe the motion of the mesh and material particles separately with respect to a reference domain.ALE has been applied successfully in many fields and particularly in metal forming processes.Despite the intensive development of these techniques, material mixing in FSW process remains very difficult to achieve numerically and especially when using finite element method [2,3].
An alternative solution for the simulation of this process is the use of meshless methods.Indeed, these techniques are developed for several decades to solve partial differential equations.They aim to avoid difficulties related to the mesh discretization and they have been intensively developed since the 90s.To our knowledge, smoothed particle hydrodynamics (SPH) was the first meshless method presented in 1977 by Lucia and Gingold and Monaghan for astrophysical problems [9,10].Another class of meshless methods called diffuse element method was proposed by Nayroles et al. [11,12].This method is based on local and moving least square approximation and a Galerkin discretization scheme.Belytschko et al. developed in 1994 element-free Galerkin method (EFG) [13].
The SPH method based on the notion of particles is best suited for large deformation problems [14] and the MLS approach [15] provides a better consistency for SPH.The first reason is that MLS approach uses monomial basis functions to ensure that a local polynomial expansion best fits the surrounding data points.The second reason is that there is less oscillation of the interpolating polynomial between particles with the MLS which is the source of the tensile instability [16].In addition, the boundary integrals were dropped in the derivation of the particle equations for SPH, so corrections are needed near boundaries [17].
The aim of the proposed work is to implement an implicit algorithm to simulate material mixing in friction stir welding.This algorithm is based on the continuous moving least square method (CMLS) [18] which combines the advantages of both approaches SPH and MLS.CMLS is formulated here in the Lagrangian framework, using the strong form of partial differential equations.This framework allows one to avoid the drawback related to numerical integration.Thermal and mechanical fields are directly attached to the points used for the domain discretization by the MLS shape functions [19].In this work, the boundary conditions are imposed with the boundary point collocation method [20], where collocation points are coincident with nodes along the boundary where the solutions are imposed.This allows one to easily simulate large deformations and material mixing around the welding tool.
The layout of this paper is as follows.In Section 2, we present the strong form of the governing equations.In Section 3, we detail the main idea of the CMLS method and we present space and time discretization procedures.In Section 4, we present the solution strategy adopted in this work.In Section 5, numerical application is proposed and a comparison with the industrial code Fluent is given to validate our model.This code is based on an Eulerian formulation and a finite volume discretization.Finally, some conclusions are presented in Section 6.

Governing Equations
Friction stir welding is a complex thermomechanical process that requires accurate knowledge of the relations existing between the main parameters of the process such as plunge depth, travel speed, rotation speed, thermomechanical properties, and tool geometry.In many contributions, the mixing zone is considered as a high viscosity incompressible fluid and the flow is obtained using computational fluid dynamics.In some other cases, the deformation is modeled using solid mechanics and numerical methods to solve the resulting nonlinear problem to compute the different variables as velocity, temperature, and so forth.In the present work, we consider a mechanical problem with a constitutive law that depends only on the equivalent strain rate.The resulting problem is then described by the conservation laws including the following.
Mass conservation: Motion equation: where  is the material density,  is the velocity vector,  is the stress tensor, and the gravity is neglected.The stress tensor , under this assumption, is given by where  is the hydrostatic pressure,  is the deviatoric stress tensor, and  is the unity tensor.To satisfy the condition of incompressibility numerically in (1), we introduce, in general, a pressure term penalized by a large factor noted .This means that the incompressibility condition ( 1) is replaced by a viscous law [21]: where ε is the strain rate tensor defined as FSW is a process that generates large strains in which elastic strains are negligible.The advantage of this hypothesis is that the material can then be modeled as viscoplastic material: where  is the material viscosity.Taking into account (3), (4), and (6), the stress tensor is given by For the viscosity, we choose a power law given by [22] where  is the temperature, , , and  are the material properties, and ε is the equivalent strain rate given by In this work, we focus our efforts on the mixing aspect that is difficult to achieve by the conventional methods such as finite element one.So, we assume that the temperature is constant and then the viscosity depends only on the equivalent strain rate ((, ε ) = ( ε )).Finally, the governing equations can be summarized as follows:

Neighboring points update
Coordinate update Problem (10) will be completed by boundary and initial conditions.

Principal of CMLS Technique. Moving least square (MLS)
technique belongs to the family of meshless methods.Since the pioneering work about diffuse elements presented by Nayroles et al. in 1992 [12], several contributions using this technique have been developed in many mathematical and mechanical fields [23,24].Continuous moving least square method (CMLS) is a modified version of MLS as shown in  [18].In this method, the approximation  ℎ of the function  is defined as a polynomial of order  but with nonconstant coefficients.The local approximation at the coordinate  ∈ Ω is given by where () is the basis function of the spatial coordinates and  is the number of the basis functions.The vector coefficients {()} are obtained by performing a weighted least square fit for the local approximation.This yields the continuous quadratic form where Ω represents the support domain of point , () is the nodal parameter of  at  = , and ( − , ℎ) is a weighting function with a compact support ℎ.The weight function plays an important role in the performance of the MLS approximation.The exponential function and spline functions are often used in practice; for example, we present a weighting function in the form of cubic spline [25]: where  = ( − )/ℎ.The vector coefficients {()} are determined in each point by minimizing the functional ().This minimization leads to the following set of linear equations for {()}: where {} is the vector that collects the nodal parameters of field function for all the points in the support domain.[()] is called the weighted moment matrix given by Solving {()} from ( 14) and substituting it into (11), the CMLS approximants can be defined as The approximation defined by ( 16) is difficult to use due to its integral form.An equivalent discrete form of this equation constructed by a particle approximation is proposed.In this method, the approximation of ( 16) at an interpolated point of coordinate , in discrete notation, leads to the following approximation: where  is the number of points in the support domain,   is the nodal parameter of  at  =  and Δ  is the volume associated with point .
The number of points , chosen in the support domain, should be sufficient to ensure that the matrix [()] in ( 15) is invertible, so as to provide the interpolation stability.The choice of  depends on the nodal distribution and the number of basis functions.In order to ensure the existence of [()] −1 and a well-conditioned [()], we usually let  ≫  [26,27].Unfortunately, there is no theoretical best value of , and it has to be determined by numerical experiments.
Recalling the form of the approximation defined in (11), so Finally, the CMLS shape functions () are defined as The CMLS shape function   () will be continuous in the entire global domain, as long as the weight functions () are chosen properly.Field points included in this support domain are used to perform the CMLS approximation for the unknown function at this point.
In the case where the interpolated point is a particle, the SPH method is used for the calculation of the particle volume [18,28] as follows: where  is the reference particle.

Algorithm for the Nearest Neighbors Search.
In meshless methods, search for neighboring points in the influence domain is the most expensive procedure in terms of time consumption.The search of neighboring points must be performed for each point in the studied domain and many times during the computation.The naive technique consists in performing a double loop over all collocation points leading to a very expensive algorithm.Several methods are proposed in the literature to improve the performance of these algorithms [17,29].In this work, we have implemented an algorithm well adapted to the studied problem.Indeed, we have proposed to construct a fixed grid containing the domain subdivided into   subdomains in the radial direction and   subdomains along the circumferential direction.The size of each subdomain is noted (, ), where  is the  subdomain length in the radial direction and  is the subdomain length in the circumferential direction (see Figure 1(a)).
The choice of  and  must depend on mesh density of collocation points in the domain and it can be considered as a user parameter.Note that the grid is fixed and each subdomain knows its neighboring subdomains.The procedure defining the neighbors of a given subdomain is executed only one time in the algorithm.For each subdomain, we have then the number of the subdomain and the list of its neighbors.The particular cases of the subdomains located at the domain edge are considered with a specific treatment.Each collocation point  is defined by a given number of the point, its coordinates (, ), and the number of of governing equations of the problem; (ii) no numerical integration is required; (iii) no mesh is used.
To apply this technique to the problem (10), the shape functions presented in Section 3.1 are used to approximate the velocity vector {} =  ⟨, V⟩, in the bidimensional case, at any point using a set of points in a local support domain of the considered point: where [  ] is the matrix of shape functions of point  and {  } is the velocity vector of point .Using ( 5) and ( 21), the strain rate vector can be obtained using the approximated velocity vector: where [  ] is the matrix of derivatives of shape functions of point .Substituting ( 21) and ( 22) in (10) and assembling, we obtain nonlinear system defined by where [] is the global mass matrix and [({})] is a global matrix dependent of {} which is the vector that collects the point velocities.3.4.Time Integration.Using an Euler implicit scheme, the time discretization of the previous problem leads to the following nonlinear system in terms of the new unknown { +1 }: where   represents the value of the unknown at time   = Δ.

Solution Strategy
Within the so-called iterative method [30], generally, a space discretization of this problem is carried out and the obtained discretized equations are solved iteratively.The nodal unknown is denoted by { +1 }.In the first iteration, one assumes in the term [({ +1 }] the value of { +1 } = {  } to start with.Therefore, the right hand side becomes known.Next one iterates until the convergence criterion is satisfied (see Figure 2).This process is controlled by a tolerance .The iterative process will continue until The solution of the problem ( 24) is the limit of the sequence ({ +1 1 }, { +1 2 }, . . ., { +1  }, . ..) which is the solution of the following recurrent algebraic systems:

Numerical Analysis and Applications
In this study, we present a numerical application to show the effectiveness of the proposed algorithm.The mixing zone which is subjected to high strain and strain rate is located in the neighborhood of the welding tool.That is why only a small and circular region is considered to model the mixing process leading to a reasonable computation time.The proposed algorithm is the first step to validate the material mixing  procedure in a bidimensional analysis.It can be extended to a thermomechanical analysis by introducing the energy equation.
In these applications, we consider a circular plate of radius   = 6 mm made of aluminum alloy 7075 with the viscoplastic material proprieties:  = 2.6910 10 N mm −2 ,  = −3.3155,and  = 0.1324 (see (8)).We perform a time integration by the implicit scheme and we choose Δ = 0.0001 s,  = 300 K, and  = 10 −6 .The welding tool is considered rigid with radius   = 1.5 mm.In the proposed work, we consider only the dwell phase and the following boundary and initial conditions:  (, , ) = 0 for  2 +  2 =  2  ∀ ≥ 0, where  = 100 rad s −1 is the rotation speed.Results of our approach are compared with those obtained using Fluent software.This code is based on an Eulerian formulation and a finite volume discretization.To validate the results of our algorithm, we choose two equivalent configurations between the two formulations update Lagrangian and Eulerian.The CMLS calculation is performed until a time  =   = 0.1 s.In the Eulerian formulation, unsteady calculation with a time   is performed.The two calculations use the same constitutive law (8).We can conclude that the distribution of points has an influence on the distribution of the components of the strain rate as shown by the results.Analysis of the results shows that the algorithm converges from the interpoint distance  = 2.5 ⋅ 10 −4 m, which is equivalent to 1780 points.In the following studies, the equivalent configuration to the interpoint distance  = 2.5 ⋅ 10 −4 m is used in different simulations because this is the optimum choice.

Influence of the Support Size.
The support size must be large enough to cover several collocation points and to avoid numerical drawback and it must be sufficiently small to avoid obtaining very smooth solution.In this simulation, we choose a fixed parameter ℎ.In the general case, this parameter should vary during the simulation to be adapted to the change of the point distribution in the domain.In this study, the variation of the interpoints distances is very small justifying the use of ℎ constant.
In this section, the goal is to determine the optimal size of the support (influence domain).For this, we chose a range of values of this size ℎ as shown in Table 1.The results show the existence of an optimum equal to ℎ = 2.5 for which the maximum relative difference, in comparison with Fluent results, is less than 5%.In the following studies, the optimum choice of support size ℎ = 2.5 is used in different simulations.

Influence of the Weighting Function.
The CMLS method based on strong formulation forces us to treat higher order derivatives; therefore, it is necessary to use higher order weighting function to ensure the numerical computation of these derivatives.We propose here to study the effect of these weighting functions on the effectiveness of the proposed algorithm.Table 2 shows that the quality of the solution obtained by the algorithm is better in using a weighting function of an order three (cubic spline) or more wherein the relative difference is less than 5%.

Validation of the Implicit Algorithm Using the Industrial
Code Fluent.In the present section, we propose a comparative and qualitative study of our algorithm with the industrial code Fluent using the optimum choices.
Figure 7 shows the geometry and the space distribution of the material points.Two colors are used to differentiate the two plates to be welded.
Figure 8 presents different configurations of material mixing at different times.One can observe the high deformation level around the welding tool.This result, which is difficult to achieve using finite element method, shows the potential of meshless techniques for modeling FSW process.
The distributions of  and V components of the velocity vector are shown in Figures 9 and 10, respectively.Figures 12  and 13 show the evolution of  and V along the horizontal and vertical sections and the comparison between CMLS and Fluent results (see Figure 11).One can observe that the maximum relative difference is less than 1% confirming the validity of the proposed algorithm.Figures 14,15,and 16 show the distribution of the strain rate components ( ε  , ε  , ε  ) with the comparison between MLS and Fluent results.A maximum relative difference of 5% is obtained.

Conclusion
In the present work, we have proposed an implicit iterative algorithm to simulate material mixing observed in friction stir welding which is very difficult to achieve using finite element method.It is based on the combination of continuous moving least square spatial discretization, an Euler implicit scheme time discretization, and an iterative method.The performance of this implicit algorithm has been tested on the viscoplastic problem.Results of the proposed algorithm are compared with those obtained by the industrial code Fluent.The test carried out on the considered example establishes that the proposed algorithm is robust and efficient.In this work, we considered sticking conditions between tool and plates to be welded.The contact conditions between shoulder and plates are the next challenges to solve before this algorithm can be applied to FSW.We limited ourselves to bidimensional geometry but the work is in progress to extend the proposed algorithm to three-dimensional one with a complex geometry of the welding tool.

Figure 12 :
Figure 12: Evolution of  component of the velocity vector along the horizontal section ( − ).

Figure 13 :
Figure 13: Evolution of V component of the velocity vector along the horizontal section ( − ).

Figure 14 :
Figure 14: Distribution of ε  component of the strain rate given by CMLS (a) and Fluent (b).

5. 1 .
Influence of the Distribution of Points.In this section, to verify the robustness of the proposed numerical algorithm, we propose to study at first some cases of simulation to test the algorithm with different distributions of points in the study domain.These tests allow us to perform a study on the sensitivity to the density of points.We have represented, in Figure3, the domain configuration for four different distributions of points, where ℎ = 4 and  is the interpoint distance.The results shown in these numerical tests shall be registered at time   = 0.1 s.Figures 4, 5, and 6 show the distribution of the strain rate components ε  , ε  , and ε  obtained with the CMLS algorithm at different interpoints distances  = 5 ⋅ 10 −4 m,  = 3 ⋅ 10 −4 m,  = 2.5 ⋅ 10 −4 m, and  = 2.25 ⋅ 10 −4 m respectively.

Figure 15 :Figure 16 :
Figure 15: Distribution of ε  component of the strain rate given by CMLS (a) and Fluent (b).

Table 1 :
Influence of the parameter ℎ on the results of the proposed algorithm.

Table 2 :
Influence of the weighting function on the results obtained by the proposed algorithm.