Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 486346, 16 pages

http://dx.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

^{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

Academic Editor: Vladimir Turetsky

Copyright © 2015 C. T. Wu et al. 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.

#### 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 [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 () rotational boundary condition; () frictional contact algorithm; () support for very high levels of deformation; () elastic-plastic or elastic-viscoplastic material models; and () 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 -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 -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 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 [39] with denoting the isotropic thermal conductivity. is the gradient operator with respect to current position and “·” denotes the divergence operator. The Taylor-Quinney [40] 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 [29]: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 [21] 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 [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 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 [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: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 [41]. Substituting (50) into (48) yieldsAccording to (45), the discrete mortar normal contact gap and tangential slip gap at slave node are given [41] 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 [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] 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 [39] 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:

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.