An Explicit Coupled Method of FEM and Meshless Particle Method for Simulating Transient Heat Transfer Process of Friction Stir Welding

Friction stir welding (FSW) is a favorable welding technology for aluminum alloys. +e FSW process involves complex heat and mass transfer. Explicit meshless particle methods are currently popular methods for simulating the process, but they require expensive computational cost. Coupling explicit finite element method (FEM) and meshless particle methods can ease the problem by making use of high efficiency of FEM and advantages of meshless particle methods. +ough many efforts have been made to couple FEM and meshless particle methods for transient dynamics problems, coupling them for transient heat transfer problems is seldom addressed. In this work, we focus on treating this problem. We developed an explicit coupled method of FEM and the meshless particle method presented in a previous work and used it to simulate the thermal process during FSW. In the method, FEM using lumped heat capacity matrix and low-order numerical integration is constructed to obtain high efficiency. A new coupling algorithm is proposed to link thermal calculations of the weak-form FEM and the strong-form meshless particle method. Forward Euler method is used for time integration to achieve an explicit algorithm. +e coupled method is used to calculate a numerical example having analytical solution. Calculated results show that it can achieve a good accuracy. +e method is employed to simulate FSWof Al 6061-T6 plates. It predicts thermal cycles in good agreement with experimental results. It shows an accuracy comparable to that of the meshless particle method while having a higher efficiency than the latter.


Introduction
Friction stir welding (FSW) is an advantageous joining process for aluminum alloys. It makes use of heat generated from friction between a spinning tool and workpieces to join the workpieces in a solid state. During the FSW process, complex heat and mass transfer is involved. Understanding the heat transfer and material flow processes is of great importance for industrial applications of FSW since they significantly affect joint quality.
Numerical modeling provides an effective way to understand such processes. Currently, there were many numerical approaches applied to FSW modeling, including finite element method (FEM) [1][2][3], finite difference method (FDM) [4][5][6], and meshless methods [7][8][9][10][11][12][13][14]. Due to flexibility in space discretization, convenience for adaptive refinement, high accuracy, and strong ability in modeling large deformations in a Lagrangian frame, meshless methods are good options for FSW modeling. Particularly, meshless particle methods like smoothed particle hydrodynamics (SPH) [15] possess special advantages for FSW modeling; for example, they can easily model material mixing. Meshless particle methods are a class of meshless methods that represent simulated objects with so-called particles. e particles act as not only interpolation points but also material points with various physical properties. Meshless particle methods are usually derived from strong form of governing equations and totally free of mesh, though a few of them, for example, SPH, may be derived from both strong and weak forms. ey become popular approaches for FSW modeling and attract increasing research interests in recent years [9][10][11][12][13][14]. A common disadvantage of meshless methods is their more expensive computational cost compared with FEM and FDM. Usually, meshless methods use more complicated shape functions and more interpolation points, need higher-order numerical integration (mainly for weakform meshless methods), and require time-consuming search of interpolation points in each time step for large deformation problems (for example, neighboring particle search in SPH). In addition, enforcement of boundary conditions in non-interpolating meshless methods is generally not as convenient as that in FEM.
To overcome the disadvantages of meshless methods, coupling them with FEM is a popular approach. e main idea is to use meshless methods only in the problem domain where they are required to provide higher accuracy, overcome element distortion, or achieve other special purposes, while to use FEM in the rest problem domain to ensure good efficiency. Currently, a number of studies have been performed to develop coupled methods of FEM and meshless methods. A majority of existing studies focused on coupling of FEM and weak-form meshless methods such as elementfree Galerkin (EFG) method [16][17][18], meshless local Petrov-Galerkin (MLPG) method [19,20], and meshless finite difference method [21,22]. Most of these studies [16,17,19,20] constructed interface elements and established special approximation functions in these elements to achieve a continuous transition between FEM and meshless regions.
is kind of coupling approach is preferable for accuracy, but it is inconvenient to use for three-dimensional problems. Xiao and Dhanasekar [18] introduced a collocation approach to couple FEM and EFG method. e approach does not require to construct interface elements, so it is relatively easy to use for both two-and three-dimensional problems. ere were also some studies dealing with coupling between FEM and meshless particle methods [23][24][25].
ese studies mainly developed explicit coupled methods for transient dynamics problems which only require coupling of mechanic calculations. e existing coupled methods of FEM and meshless particle methods only considered coupling of mechanic calculations between FEM and meshless particle method. Since FSW is a thermomechanical problem, thermal analysis is an important aspect for FSW modeling. us, to develop an effective coupled method of FEM and meshless particle method for FSW modeling, there is still a need to handle coupling of thermal calculations. is becomes the focus of the current work. ermal conduction in FSW is a transient process and three-dimensional problem. It involves moving heat sources, complex boundary conditions, and intensive changes of thermophysical properties of materials. It is required to develop accurate, robust, and simple coupling algorithms for its analysis.
In this work, we developed an explicit coupled method of FEM and the meshless particle method presented by Xiao et al. [12]. e coupled method is currently developed to simulate transient heat transfer during FSW process. It can be combined with existing coupled methods for transient dynamics problems to achieve an efficient thermomechanical simulation of FSW process. We proposed a simple and accurate coupling algorithm to link thermal calculations of FEM and the meshless particle method. e coupling algorithm is based on ghost particles and applying heat fluxes to element faces on coupling interfaces. It can be extended to couple FEM and other similar meshless particle methods.
e correctness and accuracy of the coupled method was validated by a numerical example with analytical solution. en, it was applied to simulation of temperature distributions during FSW of Al 6061-T6 plates. e rest of the paper is organized as follows. Section 2 gives the mathematical model for heat transfer during FSW process. Section 3 introduces formulations of FEM and the meshless particle method. Section 4 details the coupled method of FEM and the meshless particle method. Section 5 presents a numerical example to validate the accuracy of the coupled method. Section 6 gives the numerical model and results for FSW of Al 6061-T6 plates. e last section draws the main conclusions of this work.

Governing Equation.
FSW butt welding of two workpieces with the same size and material is considered in this work. A schematic diagram for the FSW butt welding is shown in Figure 1. A rotating tool with a shoulder and a pin penetrates the workpieces in the initial stage and then translates along the welding direction with a certain translational speed to weld the workpieces. Heat is generated due to the friction between the tool and the workpieces. e workpiece materials are softened by the heat and mixed under the stirring of the tool to form a joint. e heat transfer in the workpieces during the FSW process is a threedimensional transient problem. It is governed by the following equation: where In the above two equations, ρ is the density, c is the specific heat, T is the temperature, t is the time, λ is the thermal conductivity, x, y, and z are the coordinates, and Q is the volumetric heat source consisting of two parts, Q s and Q b , which account for heat generation from the tool pin's side and bottom, respectively.

Boundary Conditions.
Assuming the problem is symmetric about the weld joint, one of the two workpieces is only required to be modeled. e boundary conditions for the workpiece which is modeled are shown in Figure 1. e workpiece's sides except the one at the interface of the two workpieces are exposed to the surroundings. is is also true for the part of the workpiece's top surface which is out of the tool shoulder. Convective boundary conditions are prescribed for these surfaces: where n is the unit outward normal vector, h 1 is the convective heat transfer coefficient, T a is the ambient temperature, and T s is the temperature on the surfaces. For the workpiece's bottom surface, a convective boundary condition is also prescribed. Since the bottom surface contacts with a backing plate, an effective convective heat transfer coefficient h 2 is used to consider the heat loss from the bottom surface. e part of the workpiece's top surface which is contacting with the tool shoulder receives heat generated from friction between the tool shoulder and the workpiece. It is defined as a boundary with specified surface heat flux. Its boundary condition is given as where q sh is the surface heat flux resulting from the friction between the tool shoulder and the workpiece. e workpiece's side located at the interface of the two workpieces is given an insulated boundary condition, since the problem is symmetric about this plane, and temperature gradient normal to this plane is zero. e insulated boundary condition is given as

Heat Source
Model. In the present model, heat generation is characterized by volumetric heat sources Q s and Q b as well as the surface heat flux q sh . Q s and Q b denote heat generation from the side and bottom of the tool pin, respectively. q sh denotes the heat generation from the tool shoulder. Q s is calculated as Q b is calculated as q sh is calculated as In the above three equations, C f is the fraction of heat flowing into the workpiece, ) is the maximum shear stress for yielding where σ y (T) is the yield stress at temperature T, ω is the angular velocity of the tool, r p is the radius of the tool pin, and r is the distance from a point on the tool shoulder/tool pin's bottom to the symmetric axis of the tool. e detailed derivation for the expressions of Q s , Q b , and q sh can be found in the work of Xiao et al. [12].
As shown in Figure 2, Q s , Q b , and q sh are applied to different regions of the workpiece. Q s is applied to a hollow cylindrical region around the tool pin's side with a radial thickness of t s , and Q b is applied to a cylindrical region below the tool pin's bottom with an axial height of t b . q sh is applied to the part of workpiece's top surface contacting with the tool shoulder. q sh is only involved in the welding process, since the tool shoulder does not contact with the workpiece's top surface during the plunging process.
In the current modeling, we keep the modeled workpiece static and move the tool according to welding process parameters. With the movement of the tool, influence regions of heat sources Q s , Q b , and q sh will change, and they are determined in each time step according to the tool position and its geometric dimensions.

Formulation of FEM.
Galerkin weighted residual method is used to derive FEM formulation for transient heat transfer problem. e method requires the following residual equation, instead of strong form of the governing equation, to be satisfied: In equation (10), q is the specified heat flux. It equals q sh for the part of workpiece's top surface contacting with the tool shoulder and zero for the workpiece's side between the two workpieces. In equation (11), the convective heat transfer coefficient h equals h 1 for the workpiece's sides and top surface and equals h 2 for the workpiece's bottom surface.
In Galerkin weighted residual method, the weight functions W i ′ , W i ″ , and W ‴ i have the same form of shape function N i . e weight functions are chosen as Substituting the expressions of weight functions and residuals into the residual equation (8) gives According to divergence theorem, we have Substituting the above equation into equation (13) gives the Galerkin weak form for the transient heat transfer problem: To discretize the weak form in space, the global problem domain is divided into small elements. e change rate (time derivative) of temperature is approximated as  where n n is the total number of nodes, N i is the abovementioned shape function, and (zT i /zt) is the change rate of temperature at node. In the practical calculation of zT(x)/zt, it is not necessary to calculate shape functions of all n n nodes at point x, since there are only some nodes (4 nodes for tetrahedron element used in this work) whose shape functions at point x have nonzero values. ese nodes are the nodes of the element containing the point x.
Inserting equation (16) into equation (15) and calculating the integrals in equation (15) over elements, we obtain the discrete Galerkin weak form as where n e is the total number of elements, V e denotes the volume of an element, s e 2 denotes the element faces located on s 2 , s e 3 denotes the element faces located on s 3 , and n f2 and n f3 are the number of element faces located on s 2 and s 3 , respectively. e above equation can be rewritten as where C ij forms a consistent heat capacity matrix which has nonzero off-diagonal elements. Using such heat capacity matrix requires expensive computational cost to solve its inverse matrix. Meanwhile, it is easy to cause temperature oscillation. Lumped heat capacity matrix has advantages in overcoming the above problems, and it is often used to solve transient heat transfer problems [26,27]. Hence, the current work employs lumped heat capacity matrix. It is obtained by using nodes as integration points and uniformly lumping mass of an element to its nodes when C e ij is calculated. Since FEM shape functions possess the following property using the abovementioned method leads to where m e i � (m e /n e ) is the lumped mass of node i contributed from element e, m e is the mass of the element, and n e is the total number of nodes of the element. Substituting equation (24) into equation (19), we obtain where m i � n i e�1 m e i is the total lumped mass of node i, and n i is the total number of elements connected to node i. With the lumped heat capacity given in equation (25), equation (18) can be simplified as To calculate P 1 i , P 2 i , P 3 i , and P 4 i , it is required to choose a specific element type to determine values of involved integrands at any point in an element volume V e . In this work, standard linear four-node tetrahedron element is used. For this type of element, shape functions are volumetric coordinates. ey are linear functions of spatial coordinates, and their derivatives are constant throughout the element. In a tetrahedron element, temperature at a given point x is interpolated with nodal values as Mathematical Problems in Engineering and its spatial derivative in x-direction is interpolated as e spatial derivatives of temperature in y-and z-directions have the same form of equation (28). In equations (27) and (28), T 1 -T 4 are temperatures at the four nodes of the element, and N 1 -N 4 are shape functions of the four nodes. After the spatial derivatives of temperature are determined, components of heat flux at the point x, q x (x), q y (x), and q z (x), can be easily calculated with equation (1b). In FSW simulation, thermal conductivity is temperaturedependent. Hence, it is necessary to calculate the temperature at the point x with equation (27) and then obtain the current thermal conductivity before the components of heat flux are calculated.
Among P 1 i to P 4 i , only P 1 i and P 3 i are actually required to be calculated in the FSW simulation of this work, because FEM is applied to the problem domain away from the tool where heat flux boundary conditions and volumetric heat sources are absent. P 1 i involves integral over tetrahedral elements, while P 3 i involves integral over triangular element faces.
e two integrals are calculated by the symmetric numerical integration formulae proposed by Hammer et al. [28]. e general form of the numerical integration formulae is expressed as the summation of values of integrands at specified integration points multiplied by corresponding weights. e choosing of integration formulae for P 1 i and P 3 i should consider both accuracy and efficiency. In our current calculation, the integral of P 1 i is evaluated by one-point integration formula which is precise for linear function, and the integral of P 3 i is evaluated by three-point integration formula which is precise for quadratic function. Note that integrand of P 1 i is constant over an element in the case of constant thermal conductivity, so the one-point integration formula is precise for the calculation of P 1 i in this case. In our FSW simulation, thermal conductivity is expressed as a cubic function in terms of temperature. Due to the linear temperature distribution over an element, integrand of P 1 i is also a cubic function of coordinates over the element. ough in this case a higher-order (at least five-point) integration formula is theoretically required for a precise calculation of P 1 i , a low-order integration formula, the one-point integration formula, is still adopted.
is is acceptable since spatial variations of temperature and thermal conductivity are not very intensive in the problem domain using FEM. More importantly, this provides great efficiency advantage over high-order integration formula, which accords with the aim of developing coupled methods. Integrand of P 3 i is quadratic function over a triangular element face even in the case of constant convective heat transfer coefficient, so the three-point integration formula is used and provides a precise calculation of P 3 i .

Formulation of the Meshless Particle Method.
Meshless particle methods for transient heat transfer problems were investigated for a long time. Great interests were focused on the development of SPH method [29,30]. Since conventional SPH method suffers from low accuracy of particle approximation, efforts were continuously made to achieve improvements. Currently, many SPH variations or similar meshless particle methods [12,31,32] were developed. In this work, the meshless particle method presented by Xiao et al. [12] is used and briefed here. In the method, a particle approximation having first-order consistency is used to discretize the governing equation of heat transfer. For a given function f(x), the particle approximation is expressed as and the particle approximations of its derivatives are given as where n p is the number of neighboring particles of particle i, W ij is the smoothing function which is taken as the cubic spline function: where d ij is the distance between particles i and j, and l ij � (l i + l j )/2 where l denotes the smoothing length of particle. e smoothing length is taken as α h times particle spacing. Here, α h is a user-specified coefficient. 6 Mathematical Problems in Engineering In the method, equations (1a) and (1b) are satisfied at each particle. For the two equations, equation (30) is used to discretize the first-order derivatives, and equation (29) is used to approximate the source term Q. en, the following discrete form of governing equation is obtained: where In the method, boundary conditions are imposed with a penalty method. e boundary conditions are expressed in the same form as zT zn where n x , n y , and n z are direction cosines of outward normal vectors at boundaries, and for boundariesspecified surfaceheat flux, 0, for insulatedboundaries.
e boundary conditions are enforced at boundary particles during the calculation of spatial derivatives of temperature. is is achieved by introducing a penalty term into the weighed error function for solving spatial derivatives of temperature at boundary particles: where α p is the penalty factor, n b is the number of boundary conditions for particle i, and Minimizing the weighed error function defined by equation (37) gives particle approximations for spatial derivatives of temperature at boundary particles: e penalty factor α p is chosen as 10 8 whose reasonability has been shown by Xiao et al. [12].
Simulation of heat transfer during FSW involves intensive changes of temperature and moving heat sources, which easily cause stability problem for the meshless particle method. To overcome the problem, a smoothing algorithm is introduced to smooth change rate of temperature. en, a corrected change rate of temperature is obtained as where ε is the smoothing parameter which is taken as 1.
(zT ci /zt) is finally adopted to update temperatures of particles. Figure 3 shows a schematic diagram for coupling of FEM and the meshless particle method. On the coupling interface, nodes of FEM and particles of the meshless particle method are placed to coincide with each other.

Coupling Algorithm.
To consider effects of the element region on heat transfer in the adjacent particle region, ghost particles are generated at nodes next to the coupling interface, as shown in Figure 3.
To determine at which nodes ghost particles are generated, neighboring nodes for each particle on the coupling interface are searched. For each particle on the coupling interface, nodes within a spherical region centered at the particle and having a radius of r crit are detected. Ghost particles are generated at all these nodes except the ones on the coupling interface. e critical radius r crit should be sufficiently large to provide enough ghost particles for an accurate calculation of particle approximation. It may vary with particles according to their smoothing lengths or employ a global value for all particles. In the current work, a global r crit is Mathematical Problems in Engineering used since relatively uniform particles are used. e effect of r crit on accuracy of the coupling method will be discussed in the section of numerical validation. Physical quantities of ghost particles are determined as follows: smoothing lengths are set as α h times nodal spacing at the coupling interface; temperatures and material properties (density, thermal conductivity, and specific heat) are taken as those of the corresponding nodes; heat fluxes are calculated by equation (34). During the calculation of the meshless particle method, both ghost particles and real particles are included as neighboring particles for an interested particle, as long as they are in the effective region of smoothing function of the interested particle. e calculation of equation (33) is only required for real particles. rough the inclusion of ghost particles in the calculation of equation (33) for real particles, effects of the element region on heat transfer in the adjacent particle region are obtained. Consequently, temperature changes of real particles near the coupling interface are accurately calculated.
To consider effects of the particle region on heat transfer in the adjacent element region, heat fluxes are imposed on element faces on the coupling interface; namely, the coupling interface is treated as a boundary with specified heat fluxes in the FEM calculation. e imposed heat fluxes are used to take into account the heat flowing from the particle region to the element region. ey are determined based on heat fluxes of particles on the coupling interface. For example, the triangular element face △ijk is on the coupling interface. en, heat flux is imposed on this element face. e components of the imposed heat flux at a point x on the face are interpolated as where subscripts x, y, and z denote components of the heat flux, subscripts i, j, and k denote numbers of nodes of the element face, and subscripts i′, j′, and k′ denote numbers of the particles corresponding to the nodes (namely, the particles with the same positions of the nodes). e heat flux along the normal direction of the element face is calculated by q Δijk � − q Δijkx n x + q Δijky n y + q Δijkz n z , where n x , n y , and n z are the direction cosines of outward normal vector of the element face. e minus sign is introduced here to define heat flux going into the element region through the element face as positive value, so that its definition is in accordance with that in the heat flux boundary condition equation (3). Since the element face is a boundary with specified heat flux, contributions of the heat flux imposed on the element face to heat flow rates of nodes can be obtained in the same form of the term P 2e i in equation (22) as Since lumped heat capacity matrix is used, the change rate of temperature induced by the imposed heat flux can be then calculated as in equation (26) zT c e above formula is given for node i of the element face. e change rates of temperature for nodes j and k are calculated in the same form. A node may belong to multiple element faces which are located on the coupling interface, so the total coupling-induced change rate of temperature for a node on the coupling interface should be the summation of contributions from all the relevant element faces where △ j denotes the jth relevant element face of node i and n Δ is the total number of relevant element faces. Adding the coupling-induced change rate of temperature to equation (26) leads to the total change rate of temperature for a node on the coupling interface where e integral in the above equation is calculated by the three-point symmetric numerical integration formulae proposed by Hammer et al. [28]. e above coupling algorithm is simple in implementation. ough FEM and the meshless particle method employ approximations of first-order consistency in this work, the coupling algorithm is applicable when higherorder approximations are used. For higher-order cases, the coupling procedure is the same as that for the first-order case as shown above. One should only employ corresponding shape functions and use higher-order numerical integration schemes for the calculation of equation (48).

Time Integration.
Time integration for discrete governing equations of both FEM and the meshless particle method is achieved with a forward Euler method. It advances temperature in each time step with for FEM, and with for the meshless particle method, where Δt is the time step. e current work uses a fixed time step for time integration. e forward Euler method is conditionally stable, so the time step used for calculation must be less than a critical time step. e critical time step is estimated for three-dimensional problems as where β is a scale factor, l Emin is the minimum element size characterized by minimum altitude for tetrahedron element, and l Pmin is the minimum smoothing length. In FSW modeling, thermophysical properties of materials such as c and λ change with time due to temperature change, so the critical time step will also change with time. To ensure the fixed time step set in initialization stage is applicable for the whole modeling process, the maximum λ and the minimum c that materials may experience are used in the estimation of critical time step to obtain a relatively conservative value. e reference value for scale factor β is about 1. For FSW modeling, actually admissible value of β is usually smaller than the reference value and determined by testing values around the reference value.

Numerical Validation
A numerical example concerning heat transfer in a solid cube is used to examine the correctness and accuracy of the coupled method. e solid cube has initial temperature T(x, 0) � 1 K everywhere. It is exposed to the surroundings with temperature T a � 0 K. e cube cools with time due to convective heat transfer between its surfaces and the surroundings. Convective heat transfer coefficients for all surfaces are assumed to be h � 1(W/(m 2 · K)). e cube has a side length of 2a (a � 1 m); its center is at the coordinate origin. Its material properties are assumed to be ρ � 1(kg/m 3 ), c � 1(J/(kg · K)), and λ � 1(W/(m · K)).
ere is an analytical solution available for temperature distribution of the cube during the cooling process [33], which will be used to validate the accuracy of numerical results of the coupled method. e detailed expression for the analytical solution is given as sin n i a n i aμ i cos n i x e − kn 2 i t . (53) In equation (53), k � λ/ρc, and n i denotes the positive roots of the following equation: φ(y, t) and φ(z, t) can be determined in the same way of obtaining φ(x, t). In the actual calculation of the analytical solution, φ(x, t), φ(y, t), and φ(z, t) are taken as the first 50 terms in an ascending order of n.
A coupled model shown in Figure 4(a) is established to simulate the problem. In the model, a cubical region having side length of 1 m in the center of the cube is simulated with the meshless particle method, and the rest region of the cube is simulated with FEM. e meshless region is discretized with 1,728 (12 × 12 × 12) uniform particles.
e particle spacing is about 0.09 m. e FEM region is discretized with 61,070 elements, which have an approximate global size of 0.09 m. e total simulation time is set to 0.4 s.
Firstly, we choose the coefficient for smoothing length α h as 1.5 and the critical radius r crit for determining ghost particles as twice smoothing length (about 0.27 m) for a preliminary correctness test of the coupled method. In this circumstance, l Pmin and l Emin are 0.136 m and 0.022 m, respectively.
e critical time step Δt c is estimated as 4.8 × 10 −4 s with the scale factor β taken as 1. In the preliminary test, a relatively conservative value 2.5 × 10 −4 s is used for the time step Δt. Figure 5 shows a comparison between temperature time histories of two typical points obtained by the coupled method and analytical method. Coordinates for points 1 and 2 are (0, 0, 0) and (1, 1, 1), respectively. From Figure 5, it can be seen that the results of the coupled method are in good agreement with those of the analytical method. is demonstrates the correctness of the present coupled method for transient heat transfer analysis. To quantitatively evaluate the computational accuracy, two relative errors of temperature are defined. e maximum relative error is defined as and the global relative error is defined as where T n (x i , t) and T a (x i , t) are the numerical and analytical results of temperature at node/particle i, respectively, and n t is the total number of particles and nodes used in numerical Mathematical Problems in Engineering 9 calculation. Table 1 shows the calculated e s and e g at different time. As a reference, results of FEM and the meshless particle method are also provided in the table. e FEM model and meshless model are given in Figures 4(b) and 4(c), respectively. e total number of elements/particles in the coupled model, FEM model, and meshless model is 62,798, 62,774, and 64,000 (40 × 40 × 40), respectively, which are very close. e same time step 2.5 × 10 −4 s is also used for the calculations of FEM and the meshless particle method. From Table 1, it can be seen that during the whole simulation process the maximum values of e g and e s for the coupled method are only 0.0039 and 0.0578, respectively. is indicates that the coupling method can achieve reasonable accuracy. From Table 1, it can be also observed that the accuracy of the coupling method with the current parameters seems to be lower than that of FEM and the meshless particle method. However, this situation will be obviously improved with more reasonable parameters, which will be seen in the following parametric effects study. Figure 6(a) shows the relative error of temperature on the top surface (the plane of z � 1 m) at t � 0.4 s. is surface belongs to the FEM region. e maximum absolute value of relative error on this plane is about 0.0140. Figure 6(b) shows the relative error of temperature on the plane of z � 0.045 m at t � 0.4 s.
is plane crosses through both the meshless region and FEM region, so it is convenient to observe the relative error near coupling interfaces (dashed lines in the figure). It can be seen that the most significant error occurs at the four corners of coupling interfaces. A reason for this may be that particles at the corners of interfaces have more neighboring ghost particles whose positions are determined by FEM nodes, and the uniformity of neighboring particle distribution for these particles is worse than that for other particles, so approximations for them will be less accurate. Besides, it can be seen that the maximum absolute value of relative error on this plane is about 0.0231 and larger than that on the top surface. en, parametric effects on the accuracy of the coupled method are investigated to identify optimal parametric values. Effect of critical radius r crit is firstly studied. r crit is given as α r times smoothing length. Four values of α r are taken to perform the calculation of the problem with the coupled method. Table 2 shows the relative error for different critical distance. It can be seen that as α r increases, e s decreases obviously, and e g also decreases on the whole. is means that the accuracy of the coupled method is improved with the increase of α r . e accuracy of results for α r � 4 and α r � 5 is significantly improved compared with the result for α r � 2 and becomes close to that of the meshless result and FEM result shown in Table 1. Meanwhile, the results for these two cases have the same accuracy.
is is because smoothing function has only an influence domain with twice smoothing length in radius (see equation (32)), and extra ghost particles included by setting α r to 5 exert no influence on the calculation of real particles. Since computation cost increases with the increase of number of ghost particles, α r � 4 is a preferable choice to obtain both good accuracy and high efficiency and used in the following study. Figures 6(c) and 6(d) show the relative error of temperature on the top surface and the plane of z � 0.045 m, respectively, for the case of α r � 4. In this case, the maximum absolute value of relative  error on the plane of z � 0.045 m becomes less than that on the top surface, which is in contrast to the observation for the case of α r � 2. Moreover, as seen in Figure 6(d), the error at corners of coupling interfaces is small (only about 0.0052). It is significantly reduced compared with that for the case of α r � 2 (about 0.0231), and it is close to the error at corners of the plane which belong to FEM domain. is indicates that the local error induced by coupling calculation can be controlled at a low level and will not degrade the global solution in the case of α r � 4. Effect of time step ∆t on the accuracy of the coupled method is then studied with α r � 4 and α h � 1.5. Table 3 shows the relative error for different time step. It can be seen that the time step has insignificant effect on the computational accuracy. Beside the time steps given in the table, a larger time step 6.0 × 10 −4 s is also tested. It is found that the   time step leads to numerical instability. Effect of smoothing length on the accuracy of the coupled method is also studied with α r � 4 and ∆t � 5.0 × 10 −4 s. Table 4 shows the relative error for different smoothing length. e coupled method can achieve good accuracy for all the three smoothing lengths, and its accuracy seems to be slightly better for smaller smoothing length when smoothing length is in the range of 1.0 to 1.5 times particle spacing.

Numerical Model.
In this section, the coupled method is applied to simulate FSW butt welding of a pair of Al 6061-T6 plates (Song and Kovacevic, 2004). e two aluminum alloy plates have the same dimensions which are 254 mm (length) × 102 mm (width) × 12.7 mm (thickness). e tool is made of H-13 tool steel. e main geometric dimensions of the tool are 6 mm of pin radius (r P ), 12 mm of pin height (h P ), and 25 mm of shoulder radius (r sh ). e welding process parameters are 637 r/min of tool's rotational speed and 1.59 mm/s of tool's translational speed.
As mentioned in Section 2, a half model is established to simulate the FSW process by making use of the symmetry. e critical radius r crit for determining ghost particles is taken as 4 times smoothing length of real particles on the coupling interface, i.e., 10 mm. l Pmin and l Emin for the current coupled model are 0.85 mm and 1.25 mm, respectively. e critical time step Δt c is estimated as 8 × 10 −3 s with the scale factor β taken as 1. e actually admissible time step is found to be between 5 × 10 −3 s and 6 × 10 −3 s by trying different time steps around the estimated Δt c . e time step for the current calculation is finally set to 5 × 10 −3 s. e total simulation time is 200 s. For comparison, a meshless model is also established, as shown in Figure 7(b). e meshless particle method uses the same time step for calculation. e number of particles in the meshless model is 206,582. It is very close to the total number of particles and elements in the coupled model (206,773), so it is convenient to compare computational cost.
Constants of heat source model used in the current calculation are given in Table 5. e workpiece's density is 2,700 kg/m 3 . Temperature dependence of specific heat, thermal conductivity, and yield stress are considered for the workpiece material. e specific heat and thermal conductivity are defined as a cubic function of temperature. e temperature dependence of yield stress is defined in a       tabulated form. e parameters in temperature-dependent specific heat, thermal conductivity, and yield stress are taken with the same values as those in Xiao et al. [12]. e convective heat transfer coefficients h 1 and h 2 are 30 W/(m 2 ·K) and 200 W/(m 2 ·K), respectively. Figure 8 shows temperature distributions of the workpiece calculated by the coupled method. From this figure, it can be seen that the coupled method can stably reproduce the change of temperature field during the FSW process. Figure 8(a) presents the result for t � 2 s. At this time, the pin has partially penetrated into the workpiece, and the shoulder does not touch the workpiece's top surface yet. Heat is generated from part of the pin side and the bottom of the pin. us, the highest temperature zone appears around the pin. Figure 8(b) gives the result for t � 3 s. At this time, the pin has fully penetrated into the workpiece, and the shoulder has contacted with the workpiece's top surface. Both the pin and shoulder produce heat, but the heat generated from the shoulder becomes a dominant heat source. us, the top region of the workpiece within the shoulder becomes the highest temperature zone. is implies that the welding achieves a quasi-steady state. Temperature distributions obtained by the meshless particle method with the model shown in Figure 7(b) are given in Figure 9 for comparison. e results of the coupled method are in very good agreement with those of the meshless particle method. In addition, temperature distributions obtained by the meshless particle method with a finer model have been presented by Xiao et al. [12]. ey also agree well with those presented in the current work.   e maximum error between the peak temperatures obtained by calculation and experiment is only 35 K. From Figure 10, it can be seen that thermal cycles of the four points calculated by the coupled method and the meshless particle method are very close to each other. Figure 11 gives the change of maximum temperature of the workpiece. e result of the coupled method is in good agreement with the meshless particle method. e maximum temperature rises sharply in the tool plunge process and the initial stage of welding process. en, it shows a very slow increase tendency during most of the welding process. When the welding process reaches the end, the tool is drawn out of the workpiece, and heat input disappears. Consequently, the maximum temperature drops rapidly. e maximum temperature experienced in the workpiece is only 815.1 K. It is below the material melting point 855 K. Hence, the FSW process is a solid-state joining process.

Results and Discussion.
To compare computational cost of the coupled method and the meshless particle method, CPU time consumed for calculations is recorded based on the same computer configurations. e meshless particle method requires about 21.5 hours to complete the calculation, while the coupled method takes only about 9.3 hours to finish the calculation. e total number of time steps is the same for the two approaches, so the average computational cost of a time step for the coupled method is only about 43% of that for the meshless particle method. Note that search of neighboring particles in influence domain and calculation of complicate shape functions are only required to be done once at the first time step for transient heat transfer analysis. However, in transient thermal-mechanical analysis, these time-consuming operations are required in every time step. us, the coupled method will be more competitive over the meshless particle method for thermalmechanical analysis.

Conclusions
A coupled method of FEM and meshless particle method is developed for simulating heat transfer process during FSW.
e method uses a newly developed coupling algorithm to link the thermal calculations of the two methods. Computational accuracy and parametric effect of the coupled method are examined by a numerical example. Numerical results demonstrate that the coupled method can achieve a good accuracy. e study of parametric effect shows that its accuracy is improved with the increase of critical radius r crit for determining ghost particles till it reaches about 4 times smoothing length, and has no significant dependence on time step. e coupled method works well for smoothing length between 1.0 and 1.5 times particle spacing and seems to show some improvement in accuracy with the decrease of smoothing length in this range. e coupled method is successfully applied to the simulation of heat transfer during FSW of Al 6061-T6 plates. ermal cycles of typical points calculated by the coupled method are in good agreement with experimental results. Simulation of heat transfer process during the FSW is also performed with the meshless particle method. It is observed that the two approaches produce very close results, while the coupled method requires less computational cost than the meshless particle method. e current work can provide a reference for coupling thermal calculations of FEM and other strong-form meshless particle methods. Also, incorporating the present method into existing coupled methods for transient dynamics problems will provide an effective way for thermomechanical modeling of FSW where coupled methods may be more attractive.

Data Availability
e data used to support the findings of this study are included within the article.