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

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.

C. T. Wu 1 and Wei Hu 1 and Hui-Ping Wang 2 and Hongsheng Lu 3

Academic Editor:Vladimir Turetsky

1, Livermore Software Technology Corporation, 7374 Las Positas Road, Livermore, CA 94551, USA
2, General Motors, 30500 Mound Road, Warren, MI 48090, USA
3, Hengstar Technology Corporation, Shanghai 201203, China

Received 24 March 2015; Revised 21 June 2015; Accepted 23 June 2015; 16 July 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Friction stir welding (FSW) is an innovative welding method [1] 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 [2]. 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 [3]. In particular, understanding the temperature distributions in the work piece during FSW process is an essential subject [4] due to its effects on the material flow, grain size, residual stresses, and, subsequently, the strength.

Finite element modeling provides an effective approach to better understand and visualize the influence of tool profile and input weld parameters in the welding process. It is recognized that a correct numerical model of the FSW process should avoid any unnecessary assumptions. A list of requirements [3] for multiphysics FSW analysis code should include ( [figure omitted; refer to PDF] ) rotational boundary condition; ( [figure omitted; refer to PDF] ) frictional contact algorithm; ( [figure omitted; refer to PDF] ) support for very high levels of deformation; ( [figure omitted; refer to PDF] ) elastic-plastic or elastic-viscoplastic material models; and ( [figure omitted; refer to PDF] ) support for complex geometry. The main difficulty in finite element modeling of FSW process consists in dealing with high levels of deformations involving in the complex material flow due to frictional heating and stirring motion. The Eulerian description in Computational Fluid Dynamics (CFD) models [5-7] can be easily adapted to FSW process to circumvent the mesh distortion problem encountered in the Lagrangian formulation. In this scenario, the material is analyzed as a viscous incompressible non-Newtonian fluid flowing across the Eulerian mesh and interacting with the rotating tool. Since Eulerian representation of material flow does not capture free surfaces, additional algorithms [8] are required to introduce free surfaces and thus the sticking or sliding conditions. The numerical complexity in Eulerian formulation is also concerning the needs to model the welding heat source [9] in FSW process. The heat dissipation due to frictional contact necessitates the use of the analytical [10, 11] or particular thermal models [12] for the heat source calculation that considers the rotational speed of the tool and shear stress at the interface of tool and work piece. Recently, a solid mechanics-based Eulerian model [13] was developed as an alternative to describe the heat source and correspondingly the mechanical power in the Eulerian approach. Using this modeling technique, the relationship between the rotating speed and mechanical power can be well characterized. Another numerical challenge in Eulerian modeling of FSW process is the inherent complexity in advecting material histories. When the Eulerian description of FSW process is modeled using the flow formulation, difficulties arise in the followed elastic unloading and residual stress analyses [2, 3]. An improvement of mesh distortion problem in finite element modeling of FSW process can be obtained by an [figure omitted; refer to PDF] -adaptive remeshing scheme such as the Arbitrary Lagrangian Eulerian (ALE) method [14-16]. The ALE algorithm advances the mesh independently with material flow and makes it possible to take into account the movements of free surfaces while reducing mesh distortions. An 8-node trilinear element with reduced integration and hourglass control is usually utilized to avoid the volumetric locking problem in von Mises materials. Since advancing the ALE mesh does not involve additional degrees of freedom, it is found to be efficient in the steady-state analysis [17] of FSW process. Nevertheless, the ALE approach may fail in the transient analysis when the tool travels a long distance and the high-quality mesh cannot be maintained in the stirring zone. Unlike the [figure omitted; refer to PDF] -adaptive remeshing technique in ALE formulation, the global or combined rh-adaptive remeshing schemes [18] can be employed to effectively minimize the mesh distortion in the stirring zone. Buffa et al. [19] presented a thermomechanically coupled finite element model based on a combined rh-adaptive scheme for the FSW analysis. The main difficulty of three-dimensional mesh generation comes from the existence of silvers (highly distorted tetrahedral) with extreme dihedral angles. Although the Delaunay refinement scheme can remove silvers by adding some amount of Steiner point [20], it is still an open problem to determine a nontrivial lower bound on the dihedral angles of the output tetrahedral. In addition, the rigid-viscoplastic material assumption demands a postprocessing to extract the nodal temperature history for the calculation of residual stress field induced by the thermal cycle and the mechanical action applied by the tool.

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 [21-24]. The Moving Least-Squares (MLS) approximation in Element-Free Galerkin method [25] and the Reproducing Kernel (RK) approximation in Reproducing Kernel Particle [26] 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 [27] 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 [30] on the Dirichlet boundaries. Wu and Koishi [28] 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 [31] 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) [32] for the particulate rubber composites in the near-incompressible micromechanical analysis [33]. Recently, an immersed meshfree Galerkin formulation was proposed by Wu et al. [34] 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 [35] 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 [36] has demonstrated the meshfree capability in modeling severe material flow in FSW process, the fundamental difficulties associated with the tensile instability [37] and the spurious modes [38] 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 [25] 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 [figure omitted; refer to PDF] is a bounded polygon with disjointed boundary [figure omitted; refer to PDF] . The notation [figure omitted; refer to PDF] describes a Dirichlet boundary imposed by a temperature [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is the Neumann boundary prescribed by a heat flux. We also assume that [figure omitted; refer to PDF] and [figure omitted; refer to PDF] do not vary with time. The boundary [figure omitted; refer to PDF] denotes the contact surface with a thermal exchange due the conduction between the tool and work piece. The boundary [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] per unit deformed volume from the plastic deformation, the strong form of the thermal energy conservation equation reads [figure omitted; refer to PDF] subject to boundary conditions [figure omitted; refer to PDF] and initial condition at time [figure omitted; refer to PDF] [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the mass density and [figure omitted; refer to PDF] is the heat capacity. Equation (2) is known as the Fourier law [39] with [figure omitted; refer to PDF] denoting the isotropic thermal conductivity. [figure omitted; refer to PDF] is the gradient operator with respect to current position [figure omitted; refer to PDF] and " [figure omitted; refer to PDF] ·" denotes the divergence operator. The Taylor-Quinney [40] coefficient [figure omitted; refer to PDF] in (3) takes into account the fraction of heat generated by the plastic deformation energy dissipation. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are the deviatoric parts of Cauchy stress and the rate of plastic straining, respectively. [figure omitted; refer to PDF] in (4) is the temperature imposed on the Dirichlet boundary. [figure omitted; refer to PDF] in (5) is the normal heat flux imposed on the Neumann boundary. [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] are heat transfer coefficients for conduction, convection, and radiation, respectively. The notation [figure omitted; refer to PDF] defines the outward unit normal vector on [figure omitted; refer to PDF] . The second term on the right-hand side of (6) represents the rate of frictional energy dissipation where [figure omitted; refer to PDF] is the fraction of heat generated by the frictional contact. [figure omitted; refer to PDF] is the Cauchy contact traction and [figure omitted; refer to PDF] is the contact slip rate which represents the jump in velocity across the contact. [figure omitted; refer to PDF] is the ambient temperature and [figure omitted; refer to PDF] is the temperature of the tool. [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] : [figure omitted; refer to PDF]

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

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: [figure omitted; refer to PDF] subject to Dirichlet and Neumann boundary conditions [figure omitted; refer to PDF] together with unilateral contact conditions [figure omitted; refer to PDF] and Coulomb friction conditions [figure omitted; refer to PDF] and initial conditions [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the body force vector and [figure omitted; refer to PDF] is the Cauchy stress obtained from the constitutive law. The notation [figure omitted; refer to PDF] describes a Dirichlet boundary imposed by a displacement [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is the Neumann boundary prescribed by a surface traction [figure omitted; refer to PDF] with [figure omitted; refer to PDF] . The notation [figure omitted; refer to PDF] in (14) is the gap function defined in Section 3.2. [figure omitted; refer to PDF] denotes the outward unit normal to contact surface [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is the tangential component of contact traction [figure omitted; refer to PDF] . Equation (15) states that the vector of slip rate [figure omitted; refer to PDF] should be in the direction of tangential traction. [figure omitted; refer to PDF] is the temperature dependent friction coefficient. [figure omitted; refer to PDF] stands for the Euclidean norm. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] , such that for arbitrary variation [figure omitted; refer to PDF] , the following equation is satisfied: [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] employing the thermal weak form of (11) to find [figure omitted; refer to PDF] such that [figure omitted; refer to PDF] with [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] is an index set. [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] in the reference configuration [figure omitted; refer to PDF] is described by a mapping [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] is the spatial position of the material point [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] . Therefore, the material displacement can be defined by [figure omitted; refer to PDF]

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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] be the interior of the supp [figure omitted; refer to PDF] . We assume that [figure omitted; refer to PDF] is star-shaped with respect to a ball [figure omitted; refer to PDF] and there exists a constant [figure omitted; refer to PDF] such that [figure omitted; refer to PDF] Let [figure omitted; refer to PDF] denote the nodal support radius and assume the following overlapping condition [29]: [figure omitted; refer to PDF] Often, a Lagrangian meshfree shape function [figure omitted; refer to PDF] is associated with a meshfree node [figure omitted; refer to PDF] and nodes are distinct.

Similarly, we have the mechanical weak form of (17) to find [figure omitted; refer to PDF] such that [figure omitted; refer to PDF] with [figure omitted; refer to PDF] The material displacements are also approximated using the Lagrangian meshfree shape functions [21] as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is called the "generalized" displacement of node [figure omitted; refer to PDF] . In other words, the material displacement [figure omitted; refer to PDF] is considered as an interpolant of [figure omitted; refer to PDF] in a generalized sense. In general, conventional meshfree approximations are not interpolants; that is, [figure omitted; refer to PDF] . Subsequently, the temperature field is approximated using the Lagrangian meshfree shape functions and is given by [figure omitted; refer to PDF] We also have [figure omitted; refer to PDF] 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 [30] 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 [28] 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 [figure omitted; refer to PDF] of a node set [figure omitted; refer to PDF] defined by [figure omitted; refer to PDF] the GMF method is introduced to construct a convex approximation of a given (smooth) function [figure omitted; refer to PDF] in the form of (25) and (26) such that the shape function [figure omitted; refer to PDF] satisfies the following linear polynomial reproduction property: [figure omitted; refer to PDF] With the meshfree convex approximation, we can define a [figure omitted; refer to PDF] -conforming subspace for the approximation of displacement field to be [figure omitted; refer to PDF] . 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 [29].

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: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] where [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] . Note that we have used the chain rule [figure omitted; refer to PDF] for the spatial gradient operation in the above derivations. [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] and [figure omitted; refer to PDF] denote the current configuration of deformable slave contact surface and rigid master contact surface, respectively. Choosing any contact point [figure omitted; refer to PDF] , the gap function [figure omitted; refer to PDF] in (14) is defined with respect to [figure omitted; refer to PDF] as [figure omitted; refer to PDF] Using the above definition, the two-body mortar contact virtual work defined on the nonmortar slave contact surface [figure omitted; refer to PDF] can be written by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the mortar multiplier representing the Cauchy contact traction in (6) and (14). [figure omitted; refer to PDF] and [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are numbers of nodes on the slave and master surfaces, respectively and [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are finite element interpolation functions defined on the current configuration of piecewise contact segments of [figure omitted; refer to PDF] and master contact surface [figure omitted; refer to PDF] , respectively. The associated discrete nodal values are denoted by [figure omitted; refer to PDF] , [figure omitted; refer to PDF] , and [figure omitted; refer to PDF] . [figure omitted; refer to PDF] is considered as an appropriate projection from slave contact point [figure omitted; refer to PDF] to master surface. Substituting (47) into (46) leads to the following discrete form: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are referred to as mortar integrals: [figure omitted; refer to PDF] In order to impose the normal and tangential contact constraints, the mortar multiplier is decomposed by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is unit averaged normal vector at slave node [figure omitted; refer to PDF] . The definition of [figure omitted; refer to PDF] on a discretized contact surface can be found in the literature [41]. Substituting (50) into (48) yields [figure omitted; refer to PDF] According to (45), the discrete mortar normal contact gap [figure omitted; refer to PDF] and tangential slip gap [figure omitted; refer to PDF] at slave node [figure omitted; refer to PDF] are given [41] by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is unit tensor. The normal part of mortar contact traction [figure omitted; refer to PDF] is subject to discrete form of Kuhn-Tucker conditions from (14) via [figure omitted; refer to PDF] Imposing (54) by penalty regularization leads to the following contact pressure: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is normal penalty parameter.

In mortar contact algorithm, the penalty regularization [43] 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 [41] by [figure omitted; refer to PDF] Accordingly, 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 [figure omitted; refer to PDF] to [figure omitted; refer to PDF] , and trial frictional traction can be expressed by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is frictional penalty parameter. The frictional contact traction is then corrected based on the slip condition at [figure omitted; refer to PDF] to give [figure omitted; refer to PDF] which 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 [figure omitted; refer to PDF] 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 [39] to give [figure omitted; refer to PDF] It also suffices to integrate the mechanical equation (30) by the central difference integration algorithm to yield [figure omitted; refer to PDF] where the capacity matrix [figure omitted; refer to PDF] and mass matrix [figure omitted; refer to PDF] are advantageously replaced by the lumped matrices [figure omitted; refer to PDF] ad [figure omitted; refer to PDF] in thermal and mechanical equations, respectively, in the explicit analysis. Using the Lagrangian meshfree shape functions, the position vector for integration point is updated by [figure omitted; refer to PDF] where the integration point position [figure omitted; refer to PDF] is initially located at the centroid of the integration cell and moves with material flow to its current position [figure omitted; refer to PDF] . In contrast to the finite element method, the current position of integration point [figure omitted; refer to PDF] needs not reside in the integration cell. Following the notation in (26), [figure omitted; refer to PDF] is defined to be the "generalized" current position of node [figure omitted; refer to PDF] . Subsequently, the deformation gradient at [figure omitted; refer to PDF] is computed using (39) and (61) to give [figure omitted; refer to PDF] For 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: [figure omitted; refer to PDF]

One way to sidestep this numerical difficulty is to consider the reordering of the neighboring nodal information [45, 46] such that a reconstruction of the Lagrangian meshfree shape functions is allowed and the possibility of nonpositive Jacobin determinant in (63) is suppressed. Another way to evade the nonpositive Jacobin determinant problem is to adopt an adaptive procedure similar to the combined rh-adaptive remeshing or global remeshing techniques in the finite element method. When compared with the approach of reordering of the neighboring nodal information [45, 46], the adaptive method is able to refine the nodal density and generate accurate free surfaces for a better simulation in manufacturing problems. In this study, the adaptive method based on the concept of global remeshing is adopted and an anisotropic unstructured tetrahedral mesh is pursued to model the complex geometries evolved in the formation of stirring zone. A set of meshfree nodes is extracted from the unstructured tetrahedral mesh created by a mesh generator and used for the reconstruction of the Lagrangian meshfree shape functions. The set of contact boundary nodes can also be obtained from the generated tetrahedral mesh. While most research in adaptive error estimation has focused on the development of an accurate error estimator on linear analysis, formidable difficulties still remain in nonlinear analysis. For example, the well-known Z-Z error [47] developed based on superconvergence properties for linear problems cannot quantitatively be estimated in nonlinear problem. Therefore Z-Z error estimator can only be considered as an error indicator in nonlinear problems. In engineering practice the use of an efficient error indicator turns out to be more desirable particularly in the explicit time integration method. In this study, a pointwise error indicator based on the shear deformation is used to trigger the adaptive procedure for the von Mises material. The global mesh generation in this study is comprised of an anisotropic surface Delaunay triangulation [48] and an Advancing Front tetrahedral mesh generation [49]. An adaptive solution is then completed by mapping the solution variables between the old and new spatial discretization. Apparently, the element or patch-based remap procedures [50] are not suitable for the meshfree method. In this study, we introduce a second-order accurate projection operation based on the meshfree convex approximation [29] to transfer the solution valuables. If the global remesher fails to decompose the computation domain into tetrahedrons, an alternative meshfree adaptive procedure is taken to reconstruct the neighboring nodal information using the old discretization. A sketch of the two-way adaptive procedure is illustrated in Figure 1.

Figure 1: A sketch of the two-way adaptive procedure in two-dimensional case.

[figure omitted; refer to PDF]

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 [51].

(e) Perform the anisotropic surface triangulation based on a metric map with Delaunay kernel [48] and yield the initial front of triangular faces. If there are any non-Delaunay facets in the constraining surfaces, insertion of Steiner points [52] 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 [53] 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 [figure omitted; refer to PDF] 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 by [figure omitted; refer to PDF] where [figure omitted; refer to PDF] denotes the collection of integration points at time [figure omitted; refer to PDF] . 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 = [figure omitted; refer to PDF] by [figure omitted; refer to PDF] and [figure omitted; refer to PDF] . For nodal value [figure omitted; refer to PDF] its remap is defined by the following projection operation: [figure omitted; refer to PDF] From the definition in (25), we have [figure omitted; refer to PDF] or equivalently [figure omitted; refer to PDF] where [figure omitted; refer to PDF] . Substituting (67) into (65) yields [figure omitted; refer to PDF] where [figure omitted; refer to PDF] [figure omitted; refer to PDF] in (68) is the remap function which satisfies the following interpolation property: [figure omitted; refer to PDF] as well as the following linear polynomial reproduction property: [figure omitted; refer to PDF]

Since the remap algorithm preserves linearly varying nodal values, it is second-order accurate in space. In each remap procedure, [figure omitted; refer to PDF] , [figure omitted; refer to PDF] are the newly constructed meshfree shape functions evaluated on the current configuration based on the old spatial discretization [figure omitted; refer to PDF] . The tilde symbol in (65) stands for a "generalized" nodal quantity as defined in (25). If any boundary node [figure omitted; refer to PDF] does not reside in the old spatial discretization [figure omitted; refer to PDF] , an appropriate projection is performed to find the closet point on the boundary of [figure omitted; refer to PDF] 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: [figure omitted; refer to PDF]

Similarly, another consistent remap procedure is performed for the state valuable [figure omitted; refer to PDF] such as effective plastic strain, stress components, and other internal valuables sampled at the integration point per integration cell: [figure omitted; refer to PDF] where [figure omitted; refer to PDF] The matrices [figure omitted; refer to PDF] are computed on the set of integration points defined on the current configuration using the GMF method. The summation " [figure omitted; refer to PDF] " 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 [figure omitted; refer to PDF] 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 becomes [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the decomposed deformation gradient due to the reconstruction of meshfree shape functions [figure omitted; refer to PDF] and is given by [figure omitted; refer to PDF] Here, we define [figure omitted; refer to PDF] to be a position vector on a new reference configuration [figure omitted; refer to PDF] at time [figure omitted; refer to PDF] . 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.

Figure 2: FSW simulation model.

[figure omitted; refer to PDF]

Figure 3: Dimensions of the tool and work piece.

[figure omitted; refer to PDF]

Figure 4: Initial mesh of the work piece.

[figure omitted; refer to PDF]

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.

Table 1: Material parameters.

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

Table 2: Yield stresses of the work piece.

Temperature (°C)	20	100	300	550	800	1080
[figure omitted; refer to PDF] (MPa)	324	300	253	196	131	70

Figure 5: Plunge displacement and traverse speed curve of the tool.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 [54] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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.

Figure 6: Plastic deformation contour plots.

(a) Plunge stage

[figure omitted; refer to PDF]

(b) Traverse stage

[figure omitted; refer to PDF]

Figure 7: Temperature contour plots.

(a) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

Figure 8: von Mises stress contour plots.

(a) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF] s

[figure omitted; refer to PDF]

Figure 9: Work piece energy curves and tool contact force curve.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

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 [3] 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.

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.

Conflict of Interests

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

References

[1] C. J. Dawes, W. M. Thomas, "Friction stir process welds in aluminum alloys," Welding Journal , vol. 75, no. 3, pp. 41-45, 1996.

[2] O. Lorrain, J. Serri, V. Favier, H. Zahrouni, 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.

[3] D. M. Neto, 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.

[4] H. Pashazadeh, A. Masoumi, 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.

[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.

[6] P. A. Colegrove, 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.

[7] D. Jacquin, B. de-Meester, A. Simar, D. Deloison, F. Montheillet, 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.

[8] R. Nandan, G. G. Roy, 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.

[9] A. Bastier, M. H. Maitournam, F. Roger, 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.

[10] M. Song, 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.

[11] M. Z. H. Khandkar, J. A. Khan, 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.

[12] P. A. Colegrove, 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.

[13] Z. Zhang, 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.

[14] D. H. Santiago, G. Lombera, S. Urquiza, A. Cassanelli, L. A. Vedia, "Numerical modeling of welded joints by the 'Friction Stir Welding' process," Materials Research , vol. 7, no. 4, pp. 569-574, 2004.

[15] H. Schmidt, 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.

[16] M. Chiumenti, M. Cervera, C. Agelet-de-Saracibar, N. Dialami, "Numerical modeling of friction stir welding processes," Computer Methods in Applied Mechanics and Engineering , vol. 254, pp. 353-369, 2013.

[17] Y. J. Chao, X. Qi, 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.

[18] H. Askes, A. Rodriguez-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.

[19] G. Buffa, J. Hua, R. Shivpuri, 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.

[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.

[21] J.-S. Chen, C. Pan, C.-T. Wu, 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.

[22] C.-T. Wu, J.-S. Chen, L. Chi, F. Huck, "Lagrangian meshfree formulation for analysis of geotechnical materials," Journal of Engineering Mechanics , vol. 127, no. 5, pp. 440-449, 2001.

[23] S. Li, W. K. Liu Meshfree Particle Method , Springer, Berlin, Germany, 2004.

[24] D. C. Simkins, 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, 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, Y. F. Zhang, "Reproducing kernel particle methods," International Journal for Numerical Methods in Fluids , vol. 20, no. 8-9, pp. 1081-1106, 1995.

[27] L. B. Lucy, "A numerical approach to the testing of the fission hypothesis," The Astronomical Journal , vol. 82, pp. 1013-1024, 1977.

[28] C. T. Wu, 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.

[29] C. T. Wu, C. K. Park, 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.

[30] H.-P. Wang, C.-T. Wu, Y. Guo, 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.

[31] C. K. Park, C. T. Wu, 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.

[32] C. T. Wu, W. Hu, 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.

[33] C. T. Wu, 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.

[34] C. T. Wu, Y. Guo, 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.

[35] C. T. Wu, 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.

[36] W. Pan, D. Li, A. M. Tartakovsky, S. Ahzi, M. Khraisheh, 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.

[37] C. T. Dyka, P. W. Randles, 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.

[38] T. Belytschko, Y. Guo, W. K. Liu, 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.

[39] T. J. R. Hughes The Finite Element Method , Prentice Hall, Englewood Cliffs, NJ, USA, 2000.

[40] T. D. Marusich, M. Ortiz, "Modelling and simulation of high-speed machining," International Journal for Numerical Methods in Engineering , vol. 38, no. 21, pp. 3675-3694, 1995.

[41] M. A. Puso, 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.

[42] W. Hu, 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.

[43] N. Kikuchi, 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, 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, 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.

[46] C. T. Wu, 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.

[47] O. C. Zienkiewicz, 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, H. Borouchaki Delaunay Triangulation and Meshing: Application to Finite Elements , Hermes, Paris, France, 1998.

[49] R. Lohner, 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.

[50] A. Bucher, A. Meyer, U.-J. Görke, 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.

[51] G. Farin, "Triangular Bernstein-Bezier patches," Computer Aided Geometric Design , vol. 3, no. 2, pp. 83-127, 1986.

[52] P. Moller, 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.

[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.