Research Article | Open Access

Volume 2015 |Article ID 486346 | https://doi.org/10.1155/2015/486346

C. T. Wu, Wei Hu, Hui-Ping Wang, Hongsheng Lu, "A Robust Numerical Procedure for the Thermomechanical Flow Simulation of Friction Stir Welding Process Using an Adaptive Element-Free Galerkin Method", Mathematical Problems in Engineering, vol. 2015, Article ID 486346, 16 pages, 2015. https://doi.org/10.1155/2015/486346

# A Robust Numerical Procedure for the Thermomechanical Flow Simulation of Friction Stir Welding Process Using an Adaptive Element-Free Galerkin Method

Revised21 Jun 2015
Accepted23 Jun 2015
Published16 Jul 2015

#### Abstract

A meshfree modeling technique of material flow in the three-dimensional multiphysics thermomechanical friction stir welding process is presented. In this numerical model, the discretization in space is derived by the Element-Free Galerkin method using a Lagrangian meshfree convex approximation. The discrete thermal and mechanical equations are weakly coupled as the time advances using a forward difference scheme. A mortar contact algorithm is employed to model the stirring effect and heat generation due to frictional contact. Heat conductance between contacting bodies is considered as a function of contact pressure. A two-way adaptive procedure is introduced to the coupled thermomechanical system to surpass potential numerical problems associated with the extensive material deformation and spatial discretization. In each adaptive phase, a consistent projection operation utilizing the first-order meshfree convex approximation is performed to remap the solution variables. Finally, a three-dimensional multiphysics thermomechanical coupled friction stir welding problem is analyzed to demonstrate the effectiveness of the present meshfree numerical procedure.

#### 1. Introduction

Friction stir welding (FSW) is an innovative welding method  that is viable for jointing aluminum alloys, cooper, magnesium, and other low-melting-point metallic materials. It combines frictional heating and stirring motion to soften and mix the interface between two work pieces yielding a solid and fully consolidated weld. Since FSW is a solid state jointing process, a high quality weld can be achieved with the absence of solidification cracking, porosity, oxidation, and other defects typical to traditional fusion welding. FSW offers high levels of repeatability, limited energy consumption, and ease of automation and thus gains its popularity in aerospace, automotive, railway, and nuclear industries. The FSW process consists of four basic phases, namely, the plunging, stirring, welding, and retraction. The plunging phase starts with plunging the rigid cylindrical spinning tool into the joint line until the shoulder contacts the top surface of work piece. In stirring phase, the heat is generated by means of work from frictional contact and material plastic deformation. This heat is dissipated into the welding zone and results in an increase of temperature and material softening. In particular, the heat generated from frictional contact plays an important role in determining the material rheological behavior and the success of the deposition welding process. The welding phase is initiated by moving either the tool or the work piece and enables the materials of two work pieces to mix together. Finally, the welding process stops and the tool retracts from the work piece. These four FSW phases constitute complicated multiphysics thermomechanical conditions which are very hard to determine experimentally . Although many successful experimental investigations have already been conducted to the adjustment of tool profile and input weld parameters such as tool speed, feed rate, and tool depth, several aspects of the FSW are still poorly understood and require further study . In particular, understanding the temperature distributions in the work piece during FSW process is an essential subject  due to its effects on the material flow, grain size, residual stresses, and, subsequently, the strength.

An alternative approach to model large deformation problem is the use of meshfree methods. Meshfree methods based on the Galerkin approach have been shown to effectively deal with the von Mises and pressure dependent materials in large deformation analyses . The Moving Least-Squares (MLS) approximation in Element-Free Galerkin method  and the Reproducing Kernel (RK) approximation in Reproducing Kernel Particle  are two commonly used approximations in meshfree methods. The MLS approximation was achieved by performing a position dependent weighted least squares fit and the RK approximation was achieved by introducing a correction function to the Smooth Particle Hydrodynamics (SPH) approximation  to impose the reproduction of basis functions. However, both approximations do not bear the Kronecker-delta property, which is a convenient property for imposing the Dirichlet boundary conditions. A different way of constructing an effective meshfree approximation for the solid mechanic analysis is to incorporate the convex approximants [28, 29]. The meshfree convex approximation guarantees the unique solution inside a convex hull with a minimum distributed data set and admits a weak Kronecker-delta property at the boundaries and therefore avoids the special treatments  on the Dirichlet boundaries. Wu and Koishi  provided a unified approach that can generate specific convex approximations and reproduce several existing meshfree approximations which are referred to as the Generalized Meshfree (GMF) approximations. From the dispersion analysis, it has been shown  that meshfree convex approximation exhibits smaller lagging phase and amplitude errors than conventional MLS and RK approximations in the full discretization of the wave equation. The eigenvalue analysis also reveals that a larger critical time step can be used in the explicit dynamic analysis. The meshfree convex approximation has been introduced to the low-order finite element to yield an inf-sup stable Meshfree-Enriched Finite Element (ME-FEM)  for the particulate rubber composites in the near-incompressible micromechanical analysis . Recently, an immersed meshfree Galerkin formulation was proposed by Wu et al.  using the nonconforming meshfree convex approximation to solve the interface constraint problems in composites and was shown to satisfy an optimal error estimates in the energy norm. Most recently, a general framework of the meshfree-enriched multiscale finite element formulation was developed by Wu and Hu  for the acoustic wave propagation analysis. In their multiscale finite element model, the fine-scale meshfree convex approximation satisfies the homogenous Dirichlet boundary condition and requires no a priori knowledge of fundamental solution for approximating the oscillatory type of Helmholtz equations. While those meshfree methods using meshfree convex approximation provide high fidelity predictions of material responses in the solid mechanics analyses, their robustness has not yet been fully investigated in the severe deformation analysis of thermomechanical problems such as FSW process. Although an alternative solution using SPH method  has demonstrated the meshfree capability in modeling severe material flow in FSW process, the fundamental difficulties associated with the tensile instability  and the spurious modes  remain challenging numerical issues in solid applications. In particular, the explicit dynamic nature of SPH formulation and the nontrivial essential boundary condition treatments prohibit the adoption of an implicit version for the postwelding residual stress analysis.

The objective of this study is to present a robust numerical procedure based on the Element-Free Galerkin method  using the Lagrangian meshfree convex approximation and a two-way adaptive procedure for the multiphysics simulation of FSW process. The reminder of the paper is organized as follows: the governing equations in FSW problem and their variational formulations are given in Section 2. Details on the spatial discretization of the variational formulations using Lagrangian meshfree convex approximations are given in Section 3. The numerical treatment of contact condition is also discussed in the same section. Section 4 describes the time discretization for the coupled thermomechanical equations. The computational strategies for the two-way adaptive procedure and the remap algorithm are given in Section 5. Finally, a numerical example is presented in Section 6 and conclusions are made in Section 7.

#### 2. Governing Equations

##### 2.1. Thermal Problem

We consider the transient heat transfer response of a FSW work piece in a three-dimensional case. We assume that the domain of the work piece is a bounded polygon with disjointed boundary . The notation describes a Dirichlet boundary imposed by a temperature . is the Neumann boundary prescribed by a heat flux. We also assume that and do not vary with time. The boundary denotes the contact surface with a thermal exchange due the conduction between the tool and work piece. The boundary is the surface with the thermal exchange due to convection and radiation. We further assume that the heat generation in the work piece is only due to the plastic deformation and frictional contact between tool and work piece. Giving an internal heat generation rate per unit deformed volume from the plastic deformation, the strong form of the thermal energy conservation equation readssubject to boundary conditionsand initial condition at time where is the mass density and is the heat capacity. Equation (2) is known as the Fourier law  with denoting the isotropic thermal conductivity. is the gradient operator with respect to current position and “·” denotes the divergence operator. The Taylor-Quinney  coefficient in (3) takes into account the fraction of heat generated by the plastic deformation energy dissipation. and are the deviatoric parts of Cauchy stress and the rate of plastic straining, respectively. in (4) is the temperature imposed on the Dirichlet boundary. in (5) is the normal heat flux imposed on the Neumann boundary. , , and are heat transfer coefficients for conduction, convection, and radiation, respectively. The notation defines the outward unit normal vector on . The second term on the right-hand side of (6) represents the rate of frictional energy dissipation where is the fraction of heat generated by the frictional contact. is the Cauchy contact traction and is the contact slip rate which represents the jump in velocity across the contact. is the ambient temperature and is the temperature of the tool. in (8) is the position vector in the reference configuration. It is also necessary to match the initial condition with the Dirichlet boundary condition at time :

The variational formulation of the thermal energy conservation equation can be written to find the temperature field such that for arbitrary variation the following equation is satisfied:Integration by parts the second term of left-hand side in (10) leads to

##### 2.2. Mechanical Problem

The motion of the FSW work piece is governed by the equation of motion with the prescribed boundary and initial conditions. The mechanical problem is given as follows:subject to Dirichlet and Neumann boundary conditionstogether with unilateral contact conditionsand Coulomb friction conditionsand initial conditionswhere is the body force vector and is the Cauchy stress obtained from the constitutive law. The notation describes a Dirichlet boundary imposed by a displacement and is the Neumann boundary prescribed by a surface traction with . The notation in (14) is the gap function defined in Section 3.2. denotes the outward unit normal to contact surface . is the tangential component of contact traction . Equation (15) states that the vector of slip rate should be in the direction of tangential traction. is the temperature dependent friction coefficient. stands for the Euclidean norm. and are the initial displacement and velocity, respectively.

Subsequently, the variational equation for the mechanical problem in FSW process is formulated using the integration by part to find the displacement field , such that for arbitrary variation , the following equation is satisfied:

Now the FSW problem is stated by coupling the mechanical weak form in (17) with the thermal weak form in (11) and subjecting to the prescribed Dirichlet boundaries and initial conditions.

#### 3. Spatial Discretization

##### 3.1. Thermomechanical Equations

The standard Galerkin method is formulated on a finite dimensional space employing the thermal weak form of (11) to find such thatwithwhere and is an index set. are meshfree shape functions. In order to prevent the tensile instability caused by the Eulerian kernel, the Lagrangian kernel approach [21, 22] is considered in this development. In Lagrangian formulation, the motion of a material point originally located at a position vector in the reference configuration is described by a mapping , where is the spatial position of the material point at time . Therefore, the material displacement can be defined by

Consequently, with meshfree discretization, discrete points that carry the primary unknown variables are attached to the same set of material points throughout the course of deformation. Under this consideration, the node set is the set of nodes defined in the reference configuration. In practice, the set of meshfree nodes can be taken from the finite element nodes created by a finite element mesh generator. We also let be the interior of the supp . We assume that is star-shaped with respect to a ball and there exists a constant such thatLet denote the nodal support radius and assume the following overlapping condition :Often, a Lagrangian meshfree shape function is associated with a meshfree node and nodes are distinct.

Similarly, we have the mechanical weak form of (17) to find such thatwithThe material displacements are also approximated using the Lagrangian meshfree shape functions  aswhere is called the “generalized” displacement of node . In other words, the material displacement is considered as an interpolant of in a generalized sense. In general, conventional meshfree approximations are not interpolants; that is, . Subsequently, the temperature field is approximated using the Lagrangian meshfree shape functions and is given byWe also have in the approximation of semidiscrete thermal equations. For this reason, various numerical treatments using conforming or nonconforming methods can be chosen to impose the Dirichlet boundary conditions of the coupled thermomechanical problem. One may refer to  for a bibliographical study of those resolution techniques. In this study an alternative meshfree approximation that restores a weak Kronecker-delta property at the boundary, the meshfree convex approximation [28, 29], is utilized to allow the direct treatment of Dirichlet boundary conditions for FSW problem. We employ the GMF  method to obtain the meshfree convex approximation. The first-order convex GMF approximation is constructed using the inverse tangent basis function and the cubic spline window function is chosen to be the weight function in GMF method. Giving a convex hull of a node set defined bythe GMF method is introduced to construct a convex approximation of a given (smooth) function in the form of (25) and (26) such that the shape function satisfies the following linear polynomial reproduction property:With the meshfree convex approximation, we can define a -conforming subspace for the approximation of displacement field to be . Accordingly, the same approximation space definition applies to the temperature field. More detailed information about the derivation of GMF method and the corresponding mathematical properties can be found in .

Substituting (25) and (26) into (18) and (23) using meshfree convex approximation, the semidiscrete equations of the coupled thermomechanical problem can be expressed by the following algebraic equations:wherewhere denotes the material gradient operator. For the implementation purpose, all terms associated with the volume integration are evaluated on the reference configuration and the surface integration terms are computed on the current configuration. Using Nanson’s formula, the Cauchy stress that appears in the internal force term is transformed to the first Piola-Kirchhoff stress tensor . Note that we have used the chain rule for the spatial gradient operation in the above derivations. and are vectors of generalized nodal values for the interior nodes to be solved for the temperature and displacement fields, respectively. With the meshfree convex approximation, the unknown nodal vectors of temperature and displacement fields for boundary nodes are denoted by and , respectively. For convenience, the finite element mesh is taken as the integration cells [21, 22] for the domain integration. Each integration cell occupies an initial volume needed in the domain integration and the one-point integration rule is used for each integration cell in the computation. The integration cells also provide a set of boundary nodes information for the contact traction calculation as described in the next subsection.

##### 3.2. Thermal Contact Model

The contact traction in (33) and (36) is computed by imposing the contact constraints through the mortar contact algorithm, which is briefly summarized as below. The interested readers may refer to [41, 42] for the implementation details.

Let and denote the current configuration of deformable slave contact surface and rigid master contact surface, respectively. Choosing any contact point , the gap function in (14) is defined with respect to asUsing the above definition, the two-body mortar contact virtual work defined on the nonmortar slave contact surface can be written bywhere is the mortar multiplier representing the Cauchy contact traction in (6) and (14). and are the time dependent displacement fields of slave and rigid master contact surfaces, respectively. The discrete version of mortar contact virtual work is defined by introducing appropriate shape function expansions to the contact surface fields and mortar multiplier. The mortar contact interpolation is performed as follows:where and are numbers of nodes on the slave and master surfaces, respectively and and are finite element interpolation functions defined on the current configuration of piecewise contact segments of and master contact surface , respectively. The associated discrete nodal values are denoted by , , and . is considered as an appropriate projection from slave contact point to master surface. Substituting (47) into (46) leads to the following discrete form:where and are referred to as mortar integrals:In order to impose the normal and tangential contact constraints, the mortar multiplier is decomposed bywhere is unit averaged normal vector at slave node . The definition of on a discretized contact surface can be found in the literature . Substituting (50) into (48) yieldsAccording to (45), the discrete mortar normal contact gap and tangential slip gap at slave node are given  bywhere is unit tensor. The normal part of mortar contact traction is subject to discrete form of Kuhn-Tucker conditions from (14) viaImposing (54) by penalty regularization leads to the following contact pressure:where is normal penalty parameter.

In mortar contact algorithm, the penalty regularization  of augmented Lagrangian scheme is often used to impose the frictional contact constraints, which requires the computation of the incremental tangential slip gap. From (53), an incremental tangential slip gap which is frame indifferent is given  byAccordingly, a trial state-return map strategy is employed to determine the Coulomb frictional tractions. In the tangential contact direction, we first assume no slip from time to , and trial frictional traction can be expressed bywhere is frictional penalty parameter. The frictional contact traction is then corrected based on the slip condition at to givewhich can be used to describe the discrete stick/slip conditions of (15). When the tool surface is in contact with the work piece, the standard Fourier law cannot be used to fully describe the heat transfer phenomena. This is because the contact surfaces do not physically match perfectly leading to the gas trap between the contact surfaces. In general, this heat resistance decreases as the contact pressure increases. For this reason, the heat conductance is computed as a function of contact pressure in this study.

#### 4. Time Discretization

The highly coupled and nonlinear system in the FSW thermomechanical equations is difficult to solve by simultaneous time-stepping algorithm. The large and unsymmetric system in the fully coupled thermomechanical equations inevitably involves the convergent problem and is expensive to solve, particularly in the presence of large deformation, severe contact conditions, and contact-induced thermal shock. In the application of interest, explicit and staggered time-stepping schemes [40, 44] are considered in this study.

The thermal equation in (29) is marched through time using the forward difference algorithm  to giveIt also suffices to integrate the mechanical equation (30) by the central difference integration algorithm to yieldwhere the capacity matrix and mass matrix are advantageously replaced by the lumped matrices ad in thermal and mechanical equations, respectively, in the explicit analysis. Using the Lagrangian meshfree shape functions, the position vector for integration point is updated bywhere the integration point position is initially located at the centroid of the integration cell and moves with material flow to its current position . In contrast to the finite element method, the current position of integration point needs not reside in the integration cell. Following the notation in (26), is defined to be the “generalized” current position of node . Subsequently, the deformation gradient at is computed using (39) and (61) to giveFor the computational efficiency in explicit time integration method, the material derivatives of meshfree shape functions are always computed and stored at the first Lagrangian time step and we reuse them during the time stepping.

In the staggered time-stepping scheme, the thermomechanical system is partitioned into two phases, the isothermal mechanical phase and the rigid conduction phase, during each time increment. The isothermal mechanical phase assumes a constant temperature during the mechanical computation and the rigid conduction phase considers the constant heat generation in the thermal computation at fixed configuration. A numerical stability requirement limits the maximum time increment in the explicit method for each phase. The overall stable time increment is then defined as the smaller of the two phases. In practice, the mechanical equations always set the critical time step for stability due to much smaller time scale associated with the mechanical problem.

#### 5. Two-Way Adaptive Procedure and Remap Algorithm

As mentioned earlier in Introduction, the major difficulty arises from the numerical modeling of FSW process when it comes to large deformation simulation of the complex material flow in frictional heating and stirring motion. Strictly Lagrangian approach based on a fix mesh in finite element method experiences difficulty in dealing with mesh entanglement due to unconstrained material flows. Although the meshfree Galerkin method using Lagrangian kernel helps improve the mesh entanglement problem in finite element method, it is prohibitively extending the range of meshfree applicability to model severe deformation that goes beyond the Lagrangian description; that is, the discretized deformation mapping in (63) ceases to be injective:

The basic steps in an adaptive solution strategy are summarized as below.(a)Compute the initial solution.(b)Estimate the pointwise error indicator in solution and interactively trigger the mesh generator as error indicated.(c)Input the deformed mesh at the beginning of each adaptive step and use it to define the surface domain.(d)Approximate the surface domain by a generalized Farin algorithm .(e)Perform the anisotropic surface triangulation based on a metric map with Delaunay kernel  and yield the initial front of triangular faces. If there are any non-Delaunay facets in the constraining surfaces, insertion of Steiner points  will be made and surfaces’ edges and faces will be recovered by a series of swaps or flips. In this study we exclusively consider the nonintersecting facets as triangle.(f)Conduct the Advancing Front method which consists of filling the empty space defined by the initial front, creating tetrahedrons one at a time, and updating the front with created faces. A quadtree data structure  is utilized for the proximity search and insertion or deletion of the front points.(g)Reconstruct the neighboring nodal information using the old discretization if any step in (d)~(f) fails. Go to step (i).(h)Perform the projection operation, transfer the solution valuables, and update all nodal and internal valuables.(i)Initialize the solution and construct the approximations based on the new nodal distribution using the Lagrangian meshfree shape functions.(j)Resume the numerical solution procedure.

A customized pointwise error indicator introduced in step (b) is based on a shear deformation measure and is computed from the deformation gradient in (39). This error indicator is defined bywhere denotes the collection of integration points at time . The value 0.5 in (64) is made in average sense since the deformation gradient is generally not symmetric. When the pointwise error indicator exceeds an acceptable level, the adaptive procedure is triggered to create a new discretization in the deformed configuration. The implementation details of the anisotropic surface Delaunay triangulation and the Advancing Front tetrahedral mesh generation can be found elsewhere [48, 49] and therefore are omitted in this paper. After a new mesh is generated, the remap procedure is performed to transfer the solution valuables from the old spatial discretization to the new spatial discretization. We denote the variables before and after the each remap to be superscripted with “−” and “+,” respectively. Subsequently, we denote the unstructured tetrahedral meshes before and after remap at time = by and . For nodal value its remap is defined by the following projection operation:From the definition in (25), we haveor equivalentlywhere . Substituting (67) into (65) yieldswhere in (68) is the remap function which satisfies the following interpolation property:as well as the following linear polynomial reproduction property:

Since the remap algorithm preserves linearly varying nodal values, it is second-order accurate in space. In each remap procedure, , are the newly constructed meshfree shape functions evaluated on the current configuration based on the old spatial discretization . The tilde symbol in (65) stands for a “generalized” nodal quantity as defined in (25). If any boundary node does not reside in the old spatial discretization , an appropriate projection is performed to find the closet point on the boundary of for the subsequent remapping. We proceed to show that the above remap procedure is consistent in the sense that if the new discretization is identical to the old discretization, then all nodal quantities will remain unchanged after the remap:

Similarly, another consistent remap procedure is performed for the state valuable such as effective plastic strain, stress components, and other internal valuables sampled at the integration point per integration cell:where The matrices are computed on the set of integration points defined on the current configuration using the GMF method. The summation “” denotes the total number of integration points in the old discretization. Since the stress field is smoothed in the remap procedure, the pressure oscillation in the conventional displacement-based Element-Free Galerkin formulation can be effectively suppressed. On the other hand, it is known that the Advancing Front technique usually encounters difficulty when fronts merge. As a result, highly distorted tetrahedrons can occur in the generated mesh and greatly affects the accuracy of finite element solution. Since the construction of Lagrangian meshfree shape functions in this study does not rely on the finite element mesh, a distinct advantage of the proposed method in the adaptive procedure is its insensitivity to the existence of highly distorted tetrahedrons in the mesh.

A noteworthy addition to the proposed method is its flexibility in reconstructing the neighboring nodal information without remeshing. In conventional finite element analysis, when the global remesher is unable to generate the desired discretization at time under certain geometrical conditions, the adaptive procedure is aborted and causes termination in the simulation. With the current method, the second way for adaptive procedure will step in and replaces the global remeshing step. Under this scenario, a meshfree nodal reconstruction step is taken which maintains the material quantities for all nodal and integration points in the Lagrangian setting but requires a new neighbor search [45, 46]. Using the chain rule, the calculation for the deformation gradient becomeswhere is the decomposed deformation gradient due to the reconstruction of meshfree shape functions and is given byHere, we define to be a position vector on a new reference configuration at time . Since this meshfree nodal reconstruction step does not involve remeshing, the remap procedures are not needed. The updated deformation gradients together with the reconstructed shape functions and their derivatives are used to evaluate the discrete terms in Section 3.1. This flexibility in choosing different adaptive procedures provides an incentive for the use of adaptive Element-Free Galerkin method in large deformation analysis of manufacturing problems.

#### 6. Numerical Simulation

In this section, we consider a FSW problem as shown in Figure 2. The light blue part is the tool, and the dark blue one is the work piece which is made of aluminum. The green part is rigid to provide the global constraints. The detailed dimensions are shown in Figure 3. The initial mesh of work piece is plotted in Figure 4, which is going to be refined with respect to the contact surface curvature of the tool. The tool has a conical shoulder surface to shape the material pushed out by two sequential stages of the FSW process:(1)plunge stage: the rotating tool plunges into the work piece along vertical direction;(2)traverse stage: the rotating tool travels along horizontal direction.

The material model of the work piece is assumed to be temperature dependant ideal plasticity. The material parameters are reported in Table 1. The temperature dependant yield stresses of the work piece are listed in Table 2. The melting point of work material is around 1080°C. The tool has a constant rotating speed 125 rad/s. The plunge displacement curve and traverse speed curve can be found in Figure 5. The thermal mortar contact is defined between the tool and work piece with the Coulomb frictional coefficient 0.7. A constant temperature 20°C is applied to all the parts as initial condition.

 Density (kg·m−3) Thermal (isotropic) Young’s modulus (GPa) Poisson’s ratio Heat capacity (J·kg−1·°C) Thermal conductivity (W·m−1·K−1) Tool 7850 434 60 Rigid Work piece 2700 875 175 70 0.3
 Temperature (°C) 20 100 300 550 800 1080 (MPa) 324 300 253 196 131 70

The normalized support size of the meshfree GMF approximation is 1.1. Since the density of nodal distribution varies dramatically throughout the domain due to local refinement, the actual nodal support size for every node is adjusted according to its surrounding nodal distribution to improve overall computational performance. By using local adaptivity, the integration cell size varies between 1 mm and 8 mm. The present numerical procedure was implemented into the LS-DYNA code  and the analysis was done by MPP double precision with 4 Xeon E5520 cores. There were around 400 adaptive steps that gradually increased the number of integration cells from ~4,000 to ~130,000.

Figure 6 plots the deformation in two stages, where the top view is to the left and the central cross-section view is to the right. The material flow and free surface are well captured and represented by the local adaptive remeshing. The red zone with large effective plastic strain clearly shows the stirring region during the welding process. Figure 7 shows the temperature results of both work piece and tool at three different processing steps. The model is plotted by cutting through the central cross section to better provide the temperature distribution inside the parts. At the end of plunge stage at  s, the maximum temperature is around the front edge of the tool contact surface where the most heat sources from frictional contact are generated. In the traverse stage at  s, there is a “V” shape contour due to high gradient temperature distribution in the stirring zone, and the maximum temperature is close to the melting point. Figure 8 gives the von Mises stress results, where material softening can be clearly observed as the temperature increases from the plunge stage to traverse stage. The energy curve of the work piece and the tool contact force curve are plotted in Figure 9.

#### 7. Conclusions

The multiphysics thermomechanical complexity and severe material flow of friction stir welding process make analytical models incapable of capturing all the details needed for the satisfactory quantitative prediction of temperature and stress fields generated in the work piece. Via numerical modeling techniques, various finite element models of the friction stir welding process have been developed to help visualize and study the fundamental thermomechanical behavior of work piece and input parameters of the tool. However, there still exist great numerical challenges for the finite element simulation of friction stir welding process when the issues of handling severe mesh distortion, modeling contact and traction-free boundary conditions, performing accurate state variable remap, and maintaining the quality of adaptive mesh are simultaneously presented in the model. This paper attempts to provide an alternative approach using the adaptive Element-Free Galerkin method to overcome those numerical challenges. From the best of authors’ knowledge, this is the first paper that describes a three-dimensional adaptive Element-Free Galerkin method for the thermomechanical analysis.

In our approach, the coupled thermomechanical equations using a staggered explicit time integration scheme are modeled within a Lagrangian framework. This approach facilitates a direct incorporation of the Element-Free Galerkin formulation with the adaptive procedures to circumvent the severe mesh distortion problems. The concept of the two-way adaptive procedure is introduced to bypass the numerical difficulty caused by the abortion of adaptive mesh generation. The Lagrangian meshfree convex approximation plays a key role in simplifying the boundary condition enforcement, suppressing tensile instability, minimizing the adaptivity-induced discretization sensitivity, and offering an accurate and consistent projection operation in the remap procedure. In comparison to the existing numerical methods, the present method shares the advantages of Lagrangian approach in solid mechanics applications as well as resolves the mesh distortion problems in extremely severe deformation analyses. Subsequently, the residual stress analysis based on the present numerical procedure demands no special numerical treatments and could be easily solved by a spring-back solution of the implicit formulation which will be addressed in the near future.

The proposed method meets all the computational requirements suggested by D. M. Neto and P. Neto  for FSW analysis including (1) consideration of rotational boundary condition; (2) consideration frictional contact; (3) support for very high levels of deformation; (4) support for elastic-plastic or elastic-viscoplastic material models; and (5) support for complex geometry. A further comparison with the experimental data requires a remodeling of the problem with a more realistic FSW environment such as the tool profile, feeding angle, friction coefficients, and the temperature dependent material constants, and the results will be presented separately in another paper.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

The authors wish to thank Dr. John O. Hallquist of LSTC for his support to this research. The authors would also like to thank Dr. Art Shapiro for the technical discussions.

1. C. J. Dawes and W. M. Thomas, “Friction stir process welds in aluminum alloys,” Welding Journal, vol. 75, no. 3, pp. 41–45, 1996. View at: Google Scholar
2. O. Lorrain, J. Serri, V. Favier, H. Zahrouni, and M. El-Hadrouz, “A contribution to a critical review of friction stir welding numerical simulation,” Journal of Mechanics of Materials and Structures, vol. 4, no. 2, pp. 351–369, 2009. View at: Publisher Site | Google Scholar
3. D. M. Neto and P. Neto, “Numerical modeling of friction stir welding process: a literature review,” The International Journal of Advanced Manufacturing Technology, vol. 65, no. 1–4, pp. 115–126, 2013. View at: Publisher Site | Google Scholar
4. H. Pashazadeh, A. Masoumi, and J. Teimournezhad, “A 3D numerical model to investigate mechanical, thermal and material flow characteristics in friction stir welding of copper sheets,” International Journal of Automotive Engineering, vol. 3, no. 1, pp. 328–342, 2013. View at: Google Scholar
5. P. Ulysse, “Three-dimensional modeling of the friction stir-welding process,” International Journal of Machine Tools and Manufacture, vol. 42, no. 14, pp. 1549–1557, 2002. View at: Publisher Site | Google Scholar
6. P. A. Colegrove and H. R. Shercliff, “3-Dimensional CFD modelling of flow round a threaded friction stir welding tool profile,” Journal of Materials Processing Technology, vol. 169, no. 2, pp. 320–327, 2005. View at: Publisher Site | Google Scholar
7. D. Jacquin, B. de-Meester, A. Simar, D. Deloison, F. Montheillet, and C. Desrayaud, “A simple Eulerian thermo-mechanical modeling of friction stir welding,” Journal of Materials Processing Technology, vol. 211, no. 1, pp. 57–65, 2011. View at: Publisher Site | Google Scholar
8. R. Nandan, G. G. Roy, and T. Debroy, “Numerical simulation of three-dimensional heat transfer and plastic flow during friction stir welding,” Metallurgical and Materials Transactions A: Physical Metallurgy and Materials Science, vol. 37, no. 4, pp. 1247–1259, 2006. View at: Publisher Site | Google Scholar
9. A. Bastier, M. H. Maitournam, F. Roger, and K. Dang-Van, “Modelling of the residual state of friction stir welded plates,” Journal of Materials Processing Technology, vol. 200, no. 1–3, pp. 25–37, 2008. View at: Publisher Site | Google Scholar
10. M. Song and R. Kovacevic, “Thermal modeling of friction stir welding in a moving coordinate system, and its validation,” International Journal of Machine Tools and Manufacture, vol. 43, no. 6, pp. 605–615, 2003. View at: Publisher Site | Google Scholar
11. M. Z. H. Khandkar, J. A. Khan, and A. P. Reynolds, “Prediction of temperature distribution and thermal history during friction stir welding: input torque based model,” Science and Technology of Welding and Joining, vol. 8, no. 3, pp. 165–174, 2003. View at: Publisher Site | Google Scholar
12. P. A. Colegrove and H. R. Shercliff, “Experimental and numerical analysis of aluminium alloy 7075-T7351 friction stir welds,” Science and Technology of Welding and Joining, vol. 8, no. 5, pp. 360–368, 2003. View at: Publisher Site | Google Scholar
13. Z. Zhang and H. W. Zhang, “Solid mechanics-based Eulerian model of friction stir welding,” The International Journal of Advanced Manufacturing Technology, vol. 72, no. 9–12, pp. 1647–1653, 2014. View at: Publisher Site | Google Scholar
14. D. H. Santiago, G. Lombera, S. Urquiza, A. Cassanelli, and L. A. Vedia, “Numerical modeling of welded joints by the ‘Friction Stir Welding’ process,” Materials Research, vol. 7, no. 4, pp. 569–574, 2004. View at: Publisher Site | Google Scholar
15. H. Schmidt and J. Hattel, “A local model for the thermomechanical conditions in friction stir welding,” Modelling and Simulation in Materials Science and Engineering, vol. 13, no. 1, pp. 77–93, 2005. View at: Publisher Site | Google Scholar
16. M. Chiumenti, M. Cervera, C. Agelet-de-Saracibar, and N. Dialami, “Numerical modeling of friction stir welding processes,” Computer Methods in Applied Mechanics and Engineering, vol. 254, pp. 353–369, 2013. View at: Publisher Site | Google Scholar | MathSciNet
17. Y. J. Chao, X. Qi, and W. Tang, “Heat transfer in friction stir welding—experimental and numerical studies,” Journal of Manufacturing Science and Engineering, vol. 125, no. 1, pp. 138–145, 2003. View at: Publisher Site | Google Scholar
18. H. Askes and A. Rodríguez-Ferran, “A combined rh-adaptive scheme based on domain subdivision. Formulation and linear examples,” International Journal for Numerical Methods in Engineering, vol. 51, no. 3, pp. 253–273, 2001. View at: Publisher Site | Google Scholar
19. G. Buffa, J. Hua, R. Shivpuri, and L. Fratini, “A continuum based fem model for friction stir welding—model development,” Materials Science and Engineering A, vol. 419, no. 1-2, pp. 389–396, 2006. View at: Publisher Site | Google Scholar
20. H. Si, “Constrained Delaunay tetrahedral mesh generation and refinement,” Finite Elements in Analysis and Design, vol. 46, no. 1-2, pp. 33–46, 2010. View at: Publisher Site | Google Scholar | MathSciNet
21. J.-S. Chen, C. Pan, C.-T. Wu, and W. K. Liu, “Reproducing kernel particle methods for large deformation analysis of non-linear structures,” Computer Methods in Applied Mechanics and Engineering, vol. 139, no. 1–4, pp. 195–227, 1996. View at: Publisher Site | Google Scholar | MathSciNet
22. C.-T. Wu, J.-S. Chen, L. Chi, and F. Huck, “Lagrangian meshfree formulation for analysis of geotechnical materials,” Journal of Engineering Mechanics, vol. 127, no. 5, pp. 440–449, 2001. View at: Publisher Site | Google Scholar
23. S. Li and W. K. Liu, Meshfree Particle Method, Springer, Berlin, Germany, 2004. View at: MathSciNet
24. D. C. Simkins and S. Li, “Meshfree simulations of thermo-mechanical ductile fracture,” Computational Mechanics, vol. 38, no. 3, pp. 235–249, 2006.
25. T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994.
26. W. K. Liu, S. Jun, and Y. F. Zhang, “Reproducing kernel particle methods,” International Journal for Numerical Methods in Fluids, vol. 20, no. 8-9, pp. 1081–1106, 1995. View at: Publisher Site | Google Scholar | MathSciNet
27. L. B. Lucy, “A numerical approach to the testing of the fission hypothesis,” The Astronomical Journal, vol. 82, pp. 1013–1024, 1977. View at: Publisher Site | Google Scholar
28. C. T. Wu and M. Koishi, “A meshfree procedure for the microscopic analysis of particle-reinforced rubber compounds,” Interaction and Multi-Scale Mechanics, vol. 2, no. 2, pp. 147–169, 2009. View at: Google Scholar
29. C. T. Wu, C. K. Park, and J. S. Chen, “A generalized approximation for the meshfree analysis of solids,” International Journal for Numerical Methods in Engineering, vol. 85, no. 6, pp. 693–722, 2011. View at: Publisher Site | Google Scholar
30. H.-P. Wang, C.-T. Wu, Y. Guo, and M. E. Botkin, “A coupled meshfree/finite element method for automotive crashworthiness simulations,” International Journal of Impact Engineering, vol. 36, no. 10-11, pp. 1210–1222, 2009. View at: Publisher Site | Google Scholar
31. C. K. Park, C. T. Wu, and C. D. Kan, “On the analysis of dispersion property and stable time step in meshfree method using the generalized meshfree approximation,” Finite Elements in Analysis and Design, vol. 47, no. 7, pp. 683–697, 2011. View at: Publisher Site | Google Scholar | MathSciNet
32. C. T. Wu, W. Hu, and J. S. Chen, “A meshfree-enriched finite element method for compressible and near-incompressible elasticity,” International Journal for Numerical Methods in Engineering, vol. 90, no. 7, pp. 882–914, 2012. View at: Publisher Site | Google Scholar | MathSciNet
33. C. T. Wu and M. Koishi, “Three-dimensional meshfree-enriched finite element formulation for micromechanical hyperelastic modeling of particulate rubber composites,” International Journal for Numerical Methods in Engineering, vol. 91, no. 11, pp. 1137–1157, 2012. View at: Publisher Site | Google Scholar | MathSciNet
34. C. T. Wu, Y. Guo, and E. Askari, “Numerical modeling of composite solids using an immersed meshfree Galerkin method,” Composites Part B: Engineering, vol. 45, no. 1, pp. 1397–1413, 2013. View at: Publisher Site | Google Scholar
35. C. T. Wu and W. Hu, “Multi-scale finite element analysis of acoustic waves using global residual-free meshfree enrichments,” Interaction and Multi-scale Mechanics, vol. 6, no. 2, pp. 83–105, 2013. View at: Publisher Site | Google Scholar
36. W. Pan, D. Li, A. M. Tartakovsky, S. Ahzi, M. Khraisheh, and M. Khaleel, “A new smoothed particle hydrodynamics non-newtonian model for friction stir welding: process modeling and simulation of microstructure evolution in a magnesium alloy,” International Journal of Plasticity, vol. 48, pp. 189–204, 2013. View at: Publisher Site | Google Scholar
37. C. T. Dyka, P. W. Randles, and R. P. Ingel, “Stress points for tension instability in SPH,” International Journal for Numerical Methods in Engineering, vol. 40, no. 13, pp. 2325–2341, 1997. View at: Google Scholar
38. T. Belytschko, Y. Guo, W. K. Liu, and S. P. Xiao, “A unified stability analysis of meshless particle methods,” International Journal for Numerical Methods in Engineering, vol. 48, no. 9, pp. 1359–1400, 2000. View at: Google Scholar | MathSciNet
39. T. J. R. Hughes, The Finite Element Method, Prentice Hall, Englewood Cliffs, NJ, USA, 2000. View at: MathSciNet
40. T. D. Marusich and M. Ortiz, “Modelling and simulation of high-speed machining,” International Journal for Numerical Methods in Engineering, vol. 38, no. 21, pp. 3675–3694, 1995. View at: Publisher Site | Google Scholar
41. M. A. Puso and T. A. Laursen, “A mortar segment-to-segment contact method for large deformation solid mechanics,” Computer Methods in Applied Mechanics and Engineering, vol. 193, pp. 601–629, 2004. View at: Google Scholar
42. W. Hu and C. T. Wu, “Metal forming analysis using meshfree-enriched finite element method and mortar contact algorithm,” Interaction and Multi-Scale Mechanics, vol. 6, no. 2, pp. 237–255, 2013. View at: Google Scholar
43. N. Kikuchi and J. T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, Pa, USA, 1988.
44. S. Koric, L. C. Hibbeler, and B. G. Thomas, “Explicit coupled thermo-mechanical finite element model of steel solidification,” International Journal for Numerical Methods in Engineering, vol. 78, no. 1, pp. 1–31, 2009.
45. P. C. Guan, S. W. Chi, J. S. Chen, T. R. Slawson, and M. J. Roth, “Semi-Lagrangian reproducing kernel particle method for fragment-impact problems,” International Journal of Impact Engineering, vol. 38, no. 12, pp. 1033–1047, 2011. View at: Publisher Site | Google Scholar
46. C. T. Wu and B. Ren, “A stabilized non-ordinary state-based peridynamics for the nonlocal ductile material failure analysis in metal machining process,” Computer Methods in Applied Mechanics and Engineering, vol. 291, pp. 197–215, 2015. View at: Publisher Site | Google Scholar | MathSciNet
47. O. C. Zienkiewicz and J. Z. Zhu, “A simple error estimator and adaptive procedure for practical engineering analysis,” International Journal for Numerical Methods in Engineering, vol. 24, no. 2, pp. 337–357, 1987.
48. P. George and H. Borouchaki, Delaunay Triangulation and Meshing: Application to Finite Elements, Hermes, Paris, France, 1998. View at: MathSciNet
49. R. Lohner and P. Parikh, “Generation of three-dimensional unstructured grids by the advanced-front method,” International Journal for Numerical Methods in Fluids, vol. 8, no. 10, pp. 1135–1149, 1988. View at: Publisher Site | Google Scholar
50. A. Bucher, A. Meyer, U.-J. Görke, and R. Kreißig, “A comparison of mapping algorithms for hierarchical adaptive FEM in finite elasto-plasticity,” Computational Mechanics, vol. 39, no. 4, pp. 521–536, 2007. View at: Publisher Site | Google Scholar
51. G. Farin, “Triangular Bernstein-Bézier patches,” Computer Aided Geometric Design, vol. 3, no. 2, pp. 83–127, 1986. View at: Publisher Site | Google Scholar | MathSciNet
52. P. Moller and P. Hansbo, “On advancing front mesh generation in three dimensions,” International Journal for Numerical Methods in Engineering, vol. 38, no. 21, pp. 3551–3569, 1995. View at: Publisher Site | Google Scholar | MathSciNet
53. H. Samet, The Design and Analysis of Spatial Data Structures, Addison-Wesley, Reading, Pa, USA, 1990.
54. LS-DYNA, Keyword User's Manual, LS-DYNA, 2015.

#### More related articles

Article of the Year Award: Outstanding research contributions of 2020, as selected by our Chief Editors. Read the winning articles.