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

Authors: Yihua Xiao; Hecheng Wu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Friction stir welding (FSW) is a favorable welding technology for aluminum alloys. The 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. Though 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. The coupled method is used to calculate a numerical example having analytical solution. Calculated results show that it can achieve a good accuracy. The method is employed to simulate FSW of 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.


---

## Body

## 1. 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–3], finite difference method (FDM) [4–6], and meshless methods [7–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. The 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. They become popular approaches for FSW modeling and attract increasing research interests in recent years [9–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 weak-form 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. The 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 element-free Galerkin (EFG) method [16–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. This 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. The approach does not require to construct interface elements, so it is relatively easy to use for both two- and three-dimensional problems. There were also some studies dealing with coupling between FEM and meshless particle methods [23–25]. These studies mainly developed explicit coupled methods for transient dynamics problems which only require coupling of mechanic calculations.The 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. Thus, 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. This becomes the focus of the current work. Thermal 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]. The 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. The 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. The correctness and accuracy of the coupled method was validated by a numerical example with analytical solution. Then, it was applied to simulation of temperature distributions during FSW of Al 6061-T6 plates.The rest of the paper is organized as follows. Section2 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. The last section draws the main conclusions of this work.


## 2. Mathematical Model

### 2.1. 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 Figure1. 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. The workpiece materials are softened by the heat and mixed under the stirring of the tool to form a joint. The heat transfer in the workpieces during the FSW process is a three-dimensional transient problem. It is governed by the following equation:(1a)ρc∂T∂t=−∂qx∂x+∂qy∂y+∂qz∂z+Q,where(1b)qx=−λ∂T∂x,qy=−λ∂T∂y,qz=−λ∂T∂z.Figure 1
Schematic diagram for FSW butt welding of two workpieces.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, Qs and Qb, which account for heat generation from the tool pin’s side and bottom, respectively.


### 2.2. Boundary Conditions

Assuming the problem is symmetric about the weld joint, one of the two workpieces is only required to be modeled. The boundary conditions for the workpiece which is modeled are shown in Figure1. The workpiece’s sides except the one at the interface of the two workpieces are exposed to the surroundings. This 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:(2)λ∂T∂n=h1Ta−Ts,where n is the unit outward normal vector, h1 is the convective heat transfer coefficient, Ta is the ambient temperature, and Ts 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 h2 is used to consider the heat loss from the bottom surface.The 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(3)λ∂T∂n=qsh,where qsh is the surface heat flux resulting from the friction between the tool shoulder and the workpiece.The 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. The insulated boundary condition is given as(4)λ∂T∂n=0.


### 2.3. Heat Source Model

In the present model, heat generation is characterized by volumetric heat sourcesQs and Qb as well as the surface heat flux qsh. Qs and Qb denote heat generation from the side and bottom of the tool pin, respectively. qsh denotes the heat generation from the tool shoulder. Qs is calculated as(5)Qs=2CfτTωrp2rp+ts2−rp2.Qb is calculated as(6)Qb=CfτTωrtb.qsh is calculated as(7)qsh=CfτTωr.In the above three equations,Cf is the fraction of heat flowing into the workpiece, τT=σyT/3 is the maximum shear stress for yielding where σyT is the yield stress at temperature T, ω is the angular velocity of the tool, rp 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. The detailed derivation for the expressions of Qs, Qb, and qsh can be found in the work of Xiao et al. [12].As shown in Figure2, Qs, Qb, and qsh are applied to different regions of the workpiece. Qs is applied to a hollow cylindrical region around the tool pin’s side with a radial thickness of ts, and Qb is applied to a cylindrical region below the tool pin’s bottom with an axial height of tb. qsh is applied to the part of workpiece’s top surface contacting with the tool shoulder. qsh is only involved in the welding process, since the tool shoulder does not contact with the workpiece’s top surface during the plunging process.Figure 2
Schematic diagram for applying heat sources in the workpiece material: (a) plunging process; (b) welding process.
(a)(b)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 sourcesQs, Qb, and qsh will change, and they are determined in each time step according to the tool position and its geometric dimensions.


## 2.1. 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 Figure1. 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. The workpiece materials are softened by the heat and mixed under the stirring of the tool to form a joint. The heat transfer in the workpieces during the FSW process is a three-dimensional transient problem. It is governed by the following equation:(1a)ρc∂T∂t=−∂qx∂x+∂qy∂y+∂qz∂z+Q,where(1b)qx=−λ∂T∂x,qy=−λ∂T∂y,qz=−λ∂T∂z.Figure 1
Schematic diagram for FSW butt welding of two workpieces.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, Qs and Qb, which account for heat generation from the tool pin’s side and bottom, respectively.


## 2.2. Boundary Conditions

Assuming the problem is symmetric about the weld joint, one of the two workpieces is only required to be modeled. The boundary conditions for the workpiece which is modeled are shown in Figure1. The workpiece’s sides except the one at the interface of the two workpieces are exposed to the surroundings. This 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:(2)λ∂T∂n=h1Ta−Ts,where n is the unit outward normal vector, h1 is the convective heat transfer coefficient, Ta is the ambient temperature, and Ts 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 h2 is used to consider the heat loss from the bottom surface.The 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(3)λ∂T∂n=qsh,where qsh is the surface heat flux resulting from the friction between the tool shoulder and the workpiece.The 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. The insulated boundary condition is given as(4)λ∂T∂n=0.


## 2.3. Heat Source Model

In the present model, heat generation is characterized by volumetric heat sourcesQs and Qb as well as the surface heat flux qsh. Qs and Qb denote heat generation from the side and bottom of the tool pin, respectively. qsh denotes the heat generation from the tool shoulder. Qs is calculated as(5)Qs=2CfτTωrp2rp+ts2−rp2.Qb is calculated as(6)Qb=CfτTωrtb.qsh is calculated as(7)qsh=CfτTωr.In the above three equations,Cf is the fraction of heat flowing into the workpiece, τT=σyT/3 is the maximum shear stress for yielding where σyT is the yield stress at temperature T, ω is the angular velocity of the tool, rp 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. The detailed derivation for the expressions of Qs, Qb, and qsh can be found in the work of Xiao et al. [12].As shown in Figure2, Qs, Qb, and qsh are applied to different regions of the workpiece. Qs is applied to a hollow cylindrical region around the tool pin’s side with a radial thickness of ts, and Qb is applied to a cylindrical region below the tool pin’s bottom with an axial height of tb. qsh is applied to the part of workpiece’s top surface contacting with the tool shoulder. qsh is only involved in the welding process, since the tool shoulder does not contact with the workpiece’s top surface during the plunging process.Figure 2
Schematic diagram for applying heat sources in the workpiece material: (a) plunging process; (b) welding process.
(a)(b)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 sourcesQs, Qb, and qsh will change, and they are determined in each time step according to the tool position and its geometric dimensions.


## 3. Formulation of FEM and the Meshless Particle Method

### 3.1. Formulation of FEM

Galerkin weighted residual method is used to derive FEM formulation for transient heat transfer problem. The method requires the following residual equation, instead of strong form of the governing equation, to be satisfied:(8)∫VWi′RVdV+∫s2Wi″Rs2ds+∫s3Wi‴Rs3ds=0,where V denotes the problem domain, s2 denotes surfaces of the problem domain specified with heat flux boundary conditions, s3 denotes surfaces of the problem domain specified with convective boundary conditions, Wi′, Wi″, and Wi‴ are weight functions, and RV, Rs2, and Rs3 are residuals for governing equation (1), heat flux boundary conditions (equations (3) and (4)), and convective boundary condition (equation (2)), respectively,(9)RV=−∂qx∂x+∂qy∂y+∂qz∂z+Q−ρc∂T∂t,(10)Rs2=λ∂T∂n−q¯,(11)Rs3=λ∂T∂n−hTa−Ts.In equation (10), q¯ is the specified heat flux. It equals qsh 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 h1 for the workpiece’s sides and top surface and equals h2 for the workpiece’s bottom surface.In Galerkin weighted residual method, the weight functionsWi′, Wi″, and Wi‴ have the same form of shape function Ni. The weight functions are chosen as(12)Wi′=−Wi″=−Wi‴=Ni.Substituting the expressions of weight functions and residuals into the residual equation (8) gives(13)∫VNi−∂qx∂x+∂qy∂y+∂qz∂z+Q−ρc∂T∂tdV−∫s2Niλ∂T∂n−q¯ds−∫s3Niλ∂T∂n−hTa−Tsds=0.According to divergence theorem, we have(14)−∫VNi∂qx∂x+∂qy∂y+∂qz∂zdV=∫s2+s3Niλ∂T∂nds+∫Vqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV.Substituting the above equation into equation (13) gives the Galerkin weak form for the transient heat transfer problem:(15)∫VNiρc∂T∂tdV−∫Vqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV−∫VNiQdV−∫s2Niq¯ds−∫s3NihTa−Tsds=0.To discretize the weak form in space, the global problem domain is divided into small elements. The change rate (time derivative) of temperature is approximated as(16)∂Tx∂t=∑i=1nnNix∂Ti∂t.where nn is the total number of nodes, Ni is the abovementioned shape function, and ∂Ti/∂t is the change rate of temperature at node. In the practical calculation of ∂Tx/∂t, it is not necessary to calculate shape functions of all nn 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. These 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(17)∑e=1ne∫VeNiρc∑j=1nnNj∂Tj∂tdV−∑e=1ne∫Veqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV−∑e=1ne∫VeNiQdV−∑e=1nf2∫s2eNiq¯ds−∑e=1nf3∫s3eNihTa−Tsds=0,where ne is the total number of elements, Ve denotes the volume of an element, s2e denotes the element faces located on s2, s3e denotes the element faces located on s3, and nf2 and nf3 are the number of element faces located on s2 and s3, respectively. The above equation can be rewritten as(18)∑j=1nnCij∂Tj∂t=Pi1+Pi2+Pi3+Pi4,where(19)Cij=∑e=1neCije,(20)Pi1=∑e=1nePi1e,Pi2=∑e=1nf2Pi2e,Pi3=∑e=1nf3Pi3e,Pi4=∑e=1nePi4e,(21)Cije=∫VeρcNiNjdV,(22)Pi1e=∫Veqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV,Pi2e=∫s2eNiq¯ds,Pi3e=∫s3eNihTa−Tsds,Pi4e=∫VeNiQdV.Cij 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 Cije is calculated. Since FEM shape functions possess the following property(23)Nixj=0,ifi≠j,1,ifi=j,using the abovementioned method leads to(24)Cije=0,ifi≠j,mieci,ifi=j,where mie=me/ne is the lumped mass of node i contributed from element e, me is the mass of the element, and ne is the total number of nodes of the element. Substituting equation (24) into equation (19), we obtain(25)Cij=0,ifi≠j,mici,ifi=j,where mi=∑e=1nimie is the total lumped mass of node i, and ni 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(26)∂Ti∂t=1CiiPi1+Pi2+Pi3+Pi4.To calculatePi1, Pi2, Pi3, and Pi4, it is required to choose a specific element type to determine values of involved integrands at any point in an element volume Ve. In this work, standard linear four-node tetrahedron element is used. For this type of element, shape functions are volumetric coordinates. They 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(27)Tx=N1xT1+N2xT2+N3xT3+N4xT4,and its spatial derivative in x-direction is interpolated as(28)∂Tx∂x=∂N1x∂xT1+∂N2x∂xT2+∂N3x∂xT3+∂N4x∂xT4The spatial derivatives of temperature iny- and z-directions have the same form of equation (28). In equations (27) and (28), T1–T4 are temperatures at the four nodes of the element, and N1–N4 are shape functions of the four nodes. After the spatial derivatives of temperature are determined, components of heat flux at the point x, qxx, qyx, and qzx, can be easily calculated with equation (1b). In FSW simulation, thermal conductivity is temperature-dependent. 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.AmongPi1 to Pi4, only Pi1 and Pi3 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. Pi1 involves integral over tetrahedral elements, while Pi3 involves integral over triangular element faces. The two integrals are calculated by the symmetric numerical integration formulae proposed by Hammer et al. [28]. The general form of the numerical integration formulae is expressed as the summation of values of integrands at specified integration points multiplied by corresponding weights. The choosing of integration formulae for Pi1 and Pi3 should consider both accuracy and efficiency. In our current calculation, the integral of Pi1 is evaluated by one-point integration formula which is precise for linear function, and the integral of Pi3 is evaluated by three-point integration formula which is precise for quadratic function. Note that integrand of Pi1 is constant over an element in the case of constant thermal conductivity, so the one-point integration formula is precise for the calculation of Pi1 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 Pi1 is also a cubic function of coordinates over the element. Though in this case a higher-order (at least five-point) integration formula is theoretically required for a precise calculation of Pi1, a low-order integration formula, the one-point integration formula, is still adopted. This 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 Pi3 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 Pi3.


### 3.2. 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 fx, the particle approximation is expressed as(29)fi=∑j=1npfjW˜ij,and the particle approximations of its derivatives are given as(30)∂f∂xi=∑j=1npfji∂W˜ij∂x,∂f∂yi=∑j=1npfji∂W˜ij∂y,∂f∂zi=∑j=1npfji∂W˜ij∂z,where(31)W˜ij∂W˜ij∂x∂W˜ij∂y∂W˜ij∂zT=WijA−1Pj,A=∑j=1npWijPjPjT,Pj=1xjiyjizjiT,where np is the number of neighboring particles of particle i, fji=fj−fi, xji=xj−xi, yji=yj−yi, zji=zj−zi, and Wij is the smoothing function which is taken as the cubic spline function:(32)Wij=32πl¯ij3×23−dijl¯ij2+dij/l¯ij32,0≤dijl¯ij<1,2−dij/l¯ij36,1≤dijl¯ij<2,0,dijl¯ij≥2,where dij is the distance between particles i and j, and l¯ij=li+lj/2 where l denotes the smoothing length of particle. The smoothing length is taken as αh times particle spacing. Here, αh is a user-specified coefficient.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. Then, the following discrete form of governing equation is obtained:(33)∂Ti∂t=1ρici−∑j=1npqjix∂W˜ij∂x+∑j=1npqjiy∂W˜ij∂y+∑j=1npqjiz∂W˜ij∂z+∑j=1npQjW˜ij,where(34)qix=−λ∑j=1npTji∂W˜ij∂x,qiy=−λ∑j=1npTji∂W˜ij∂y,qiz=−λ∑j=1npTji∂W˜ij∂z.In the method, boundary conditions are imposed with a penalty method. The boundary conditions are expressed in the same form as(35)∂T∂n=∂T∂xnx+∂T∂yny+∂T∂znz=ϕ,where nx, ny, and nz are direction cosines of outward normal vectors at boundaries, and(36)ϕ=hTa−Tsλ,for convective boundaries,qshλ,for boundaries specified surface heat flux,0,for insulated boundaries.The boundary conditions are enforced at boundary particles during the calculation of spatial derivatives of temperature. This is achieved by introducing a penalty term into the weighed error function for solving spatial derivatives of temperature at boundary particles:(37)e=∑j=1npWijT^j−Tj2+αpli∑j=1nbnxj∂T∂xi+nyj∂T∂yi+nzj∂T∂zi−ϕj2,where αp is the penalty factor, nb is the number of boundary conditions for particle i, and(38)T^j=Ti+∂T∂xixji+∂T∂yiyji+∂T∂zizji.Minimizing the weighed error function defined by equation (37) gives particle approximations for spatial derivatives of temperature at boundary particles:(39)∂T∂xi=B21F1+B22F2+B23F3+B24F4,∂T∂yi=B31F1+B32F2+B33F3+B34F4,∂T∂zi=B41F1+B42F2+B43F3+B44F4,where(40)B=∑jnpWijPjPjT+αpli∑jnbQjQjT,F=∑jnpWijTjPj+αpli∑jnbϕjQj,Qj=0nxjnyjnzjT.The penalty factorαp is chosen as 108 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. Then, a corrected change rate of temperature is obtained as(41)∂Tci∂t=∂Ti∂t+ε∑j=1np∂Tj∂t−∂Ti∂tW˜ij,where ε is the smoothing parameter which is taken as 1. ∂Tci/∂t is finally adopted to update temperatures of particles.


## 3.1. Formulation of FEM

Galerkin weighted residual method is used to derive FEM formulation for transient heat transfer problem. The method requires the following residual equation, instead of strong form of the governing equation, to be satisfied:(8)∫VWi′RVdV+∫s2Wi″Rs2ds+∫s3Wi‴Rs3ds=0,where V denotes the problem domain, s2 denotes surfaces of the problem domain specified with heat flux boundary conditions, s3 denotes surfaces of the problem domain specified with convective boundary conditions, Wi′, Wi″, and Wi‴ are weight functions, and RV, Rs2, and Rs3 are residuals for governing equation (1), heat flux boundary conditions (equations (3) and (4)), and convective boundary condition (equation (2)), respectively,(9)RV=−∂qx∂x+∂qy∂y+∂qz∂z+Q−ρc∂T∂t,(10)Rs2=λ∂T∂n−q¯,(11)Rs3=λ∂T∂n−hTa−Ts.In equation (10), q¯ is the specified heat flux. It equals qsh 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 h1 for the workpiece’s sides and top surface and equals h2 for the workpiece’s bottom surface.In Galerkin weighted residual method, the weight functionsWi′, Wi″, and Wi‴ have the same form of shape function Ni. The weight functions are chosen as(12)Wi′=−Wi″=−Wi‴=Ni.Substituting the expressions of weight functions and residuals into the residual equation (8) gives(13)∫VNi−∂qx∂x+∂qy∂y+∂qz∂z+Q−ρc∂T∂tdV−∫s2Niλ∂T∂n−q¯ds−∫s3Niλ∂T∂n−hTa−Tsds=0.According to divergence theorem, we have(14)−∫VNi∂qx∂x+∂qy∂y+∂qz∂zdV=∫s2+s3Niλ∂T∂nds+∫Vqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV.Substituting the above equation into equation (13) gives the Galerkin weak form for the transient heat transfer problem:(15)∫VNiρc∂T∂tdV−∫Vqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV−∫VNiQdV−∫s2Niq¯ds−∫s3NihTa−Tsds=0.To discretize the weak form in space, the global problem domain is divided into small elements. The change rate (time derivative) of temperature is approximated as(16)∂Tx∂t=∑i=1nnNix∂Ti∂t.where nn is the total number of nodes, Ni is the abovementioned shape function, and ∂Ti/∂t is the change rate of temperature at node. In the practical calculation of ∂Tx/∂t, it is not necessary to calculate shape functions of all nn 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. These 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(17)∑e=1ne∫VeNiρc∑j=1nnNj∂Tj∂tdV−∑e=1ne∫Veqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV−∑e=1ne∫VeNiQdV−∑e=1nf2∫s2eNiq¯ds−∑e=1nf3∫s3eNihTa−Tsds=0,where ne is the total number of elements, Ve denotes the volume of an element, s2e denotes the element faces located on s2, s3e denotes the element faces located on s3, and nf2 and nf3 are the number of element faces located on s2 and s3, respectively. The above equation can be rewritten as(18)∑j=1nnCij∂Tj∂t=Pi1+Pi2+Pi3+Pi4,where(19)Cij=∑e=1neCije,(20)Pi1=∑e=1nePi1e,Pi2=∑e=1nf2Pi2e,Pi3=∑e=1nf3Pi3e,Pi4=∑e=1nePi4e,(21)Cije=∫VeρcNiNjdV,(22)Pi1e=∫Veqx∂Ni∂x+qy∂Ni∂y+qz∂Ni∂zdV,Pi2e=∫s2eNiq¯ds,Pi3e=∫s3eNihTa−Tsds,Pi4e=∫VeNiQdV.Cij 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 Cije is calculated. Since FEM shape functions possess the following property(23)Nixj=0,ifi≠j,1,ifi=j,using the abovementioned method leads to(24)Cije=0,ifi≠j,mieci,ifi=j,where mie=me/ne is the lumped mass of node i contributed from element e, me is the mass of the element, and ne is the total number of nodes of the element. Substituting equation (24) into equation (19), we obtain(25)Cij=0,ifi≠j,mici,ifi=j,where mi=∑e=1nimie is the total lumped mass of node i, and ni 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(26)∂Ti∂t=1CiiPi1+Pi2+Pi3+Pi4.To calculatePi1, Pi2, Pi3, and Pi4, it is required to choose a specific element type to determine values of involved integrands at any point in an element volume Ve. In this work, standard linear four-node tetrahedron element is used. For this type of element, shape functions are volumetric coordinates. They 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(27)Tx=N1xT1+N2xT2+N3xT3+N4xT4,and its spatial derivative in x-direction is interpolated as(28)∂Tx∂x=∂N1x∂xT1+∂N2x∂xT2+∂N3x∂xT3+∂N4x∂xT4The spatial derivatives of temperature iny- and z-directions have the same form of equation (28). In equations (27) and (28), T1–T4 are temperatures at the four nodes of the element, and N1–N4 are shape functions of the four nodes. After the spatial derivatives of temperature are determined, components of heat flux at the point x, qxx, qyx, and qzx, can be easily calculated with equation (1b). In FSW simulation, thermal conductivity is temperature-dependent. 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.AmongPi1 to Pi4, only Pi1 and Pi3 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. Pi1 involves integral over tetrahedral elements, while Pi3 involves integral over triangular element faces. The two integrals are calculated by the symmetric numerical integration formulae proposed by Hammer et al. [28]. The general form of the numerical integration formulae is expressed as the summation of values of integrands at specified integration points multiplied by corresponding weights. The choosing of integration formulae for Pi1 and Pi3 should consider both accuracy and efficiency. In our current calculation, the integral of Pi1 is evaluated by one-point integration formula which is precise for linear function, and the integral of Pi3 is evaluated by three-point integration formula which is precise for quadratic function. Note that integrand of Pi1 is constant over an element in the case of constant thermal conductivity, so the one-point integration formula is precise for the calculation of Pi1 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 Pi1 is also a cubic function of coordinates over the element. Though in this case a higher-order (at least five-point) integration formula is theoretically required for a precise calculation of Pi1, a low-order integration formula, the one-point integration formula, is still adopted. This 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 Pi3 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 Pi3.


## 3.2. 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 fx, the particle approximation is expressed as(29)fi=∑j=1npfjW˜ij,and the particle approximations of its derivatives are given as(30)∂f∂xi=∑j=1npfji∂W˜ij∂x,∂f∂yi=∑j=1npfji∂W˜ij∂y,∂f∂zi=∑j=1npfji∂W˜ij∂z,where(31)W˜ij∂W˜ij∂x∂W˜ij∂y∂W˜ij∂zT=WijA−1Pj,A=∑j=1npWijPjPjT,Pj=1xjiyjizjiT,where np is the number of neighboring particles of particle i, fji=fj−fi, xji=xj−xi, yji=yj−yi, zji=zj−zi, and Wij is the smoothing function which is taken as the cubic spline function:(32)Wij=32πl¯ij3×23−dijl¯ij2+dij/l¯ij32,0≤dijl¯ij<1,2−dij/l¯ij36,1≤dijl¯ij<2,0,dijl¯ij≥2,where dij is the distance between particles i and j, and l¯ij=li+lj/2 where l denotes the smoothing length of particle. The smoothing length is taken as αh times particle spacing. Here, αh is a user-specified coefficient.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. Then, the following discrete form of governing equation is obtained:(33)∂Ti∂t=1ρici−∑j=1npqjix∂W˜ij∂x+∑j=1npqjiy∂W˜ij∂y+∑j=1npqjiz∂W˜ij∂z+∑j=1npQjW˜ij,where(34)qix=−λ∑j=1npTji∂W˜ij∂x,qiy=−λ∑j=1npTji∂W˜ij∂y,qiz=−λ∑j=1npTji∂W˜ij∂z.In the method, boundary conditions are imposed with a penalty method. The boundary conditions are expressed in the same form as(35)∂T∂n=∂T∂xnx+∂T∂yny+∂T∂znz=ϕ,where nx, ny, and nz are direction cosines of outward normal vectors at boundaries, and(36)ϕ=hTa−Tsλ,for convective boundaries,qshλ,for boundaries specified surface heat flux,0,for insulated boundaries.The boundary conditions are enforced at boundary particles during the calculation of spatial derivatives of temperature. This is achieved by introducing a penalty term into the weighed error function for solving spatial derivatives of temperature at boundary particles:(37)e=∑j=1npWijT^j−Tj2+αpli∑j=1nbnxj∂T∂xi+nyj∂T∂yi+nzj∂T∂zi−ϕj2,where αp is the penalty factor, nb is the number of boundary conditions for particle i, and(38)T^j=Ti+∂T∂xixji+∂T∂yiyji+∂T∂zizji.Minimizing the weighed error function defined by equation (37) gives particle approximations for spatial derivatives of temperature at boundary particles:(39)∂T∂xi=B21F1+B22F2+B23F3+B24F4,∂T∂yi=B31F1+B32F2+B33F3+B34F4,∂T∂zi=B41F1+B42F2+B43F3+B44F4,where(40)B=∑jnpWijPjPjT+αpli∑jnbQjQjT,F=∑jnpWijTjPj+αpli∑jnbϕjQj,Qj=0nxjnyjnzjT.The penalty factorαp is chosen as 108 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. Then, a corrected change rate of temperature is obtained as(41)∂Tci∂t=∂Ti∂t+ε∑j=1np∂Tj∂t−∂Ti∂tW˜ij,where ε is the smoothing parameter which is taken as 1. ∂Tci/∂t is finally adopted to update temperatures of particles.


## 4. Coupled Method of FEM and the Meshless Particle Method

### 4.1. Coupling Algorithm

Figure3 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.Figure 3
Schematic diagram for the coupling of FEM and the meshless particle method.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 Figure3. 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 rcrit are detected. Ghost particles are generated at all these nodes except the ones on the coupling interface. The critical radius rcrit 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 rcrit is used since relatively uniform particles are used. The effect of rcrit 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. The calculation of equation (33) is only required for real particles. Through 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. The imposed heat fluxes are used to take into account the heat flowing from the particle region to the element region. They are determined based on heat fluxes of particles on the coupling interface. For example, the triangular element face △ijk is on the coupling interface. Then, heat flux is imposed on this element face. The components of the imposed heat flux at a point x on the face are interpolated as(42)qΔijkxx=Nixqi′x+Njxqj′x+Nkxqk′x,qΔijkyx=Nixqi′y+Njxqj′y+Nkxqk′y,qΔijkzx=Nixqi′z+Njxqj′z+Nkxqk′z,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). The heat flux along the normal direction of the element face is calculated by(43)qΔijk=−qΔijkxnx+qΔijkyny+qΔijkznz,where nx, ny, and nz are the direction cosines of outward normal vector of the element face. The 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 Pi2e in equation (22) as(44)PΔijkic=∫ΔijkNiqΔijkdsSince 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)(45)∂TΔijkic∂t=PΔijkicCii=∫ΔijkNiqΔijkdsCii.The above formula is given for nodei of the element face. The 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(46)∂Tic∂t=∑j=1nΔ∫ΔjNiqΔjdsCii,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(47)∂Ti∂t=1CiiPi1+Pi2+Pi3+Pi4+Pic,where(48)Pic=∑j=1nΔ∫ΔjNiqΔjds.The integral in the above equation is calculated by the three-point symmetric numerical integration formulae proposed by Hammer et al. [28].The above coupling algorithm is simple in implementation. Though FEM and the meshless particle method employ approximations of first-order consistency in this work, the coupling algorithm is applicable when higher-order 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).


### 4.2. 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(49)Tit+Δt=Tit+∂Tit∂tΔt,for FEM, and with(50)Tit+Δt=Tit+∂Tcit∂tΔt,for the meshless particle method, where Δt is the time step. The current work uses a fixed time step for time integration.The forward Euler method is conditionally stable, so the time step used for calculation must be less than a critical time step. The critical time step is estimated for three-dimensional problems as(51)Δtc=βmincρlEmin2λ,cρlPmin2λ,where β is a scale factor, lEmin is the minimum element size characterized by minimum altitude for tetrahedron element, and lPmin 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. The 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.


## 4.1. Coupling Algorithm

Figure3 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.Figure 3
Schematic diagram for the coupling of FEM and the meshless particle method.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 Figure3. 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 rcrit are detected. Ghost particles are generated at all these nodes except the ones on the coupling interface. The critical radius rcrit 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 rcrit is used since relatively uniform particles are used. The effect of rcrit 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. The calculation of equation (33) is only required for real particles. Through 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. The imposed heat fluxes are used to take into account the heat flowing from the particle region to the element region. They are determined based on heat fluxes of particles on the coupling interface. For example, the triangular element face △ijk is on the coupling interface. Then, heat flux is imposed on this element face. The components of the imposed heat flux at a point x on the face are interpolated as(42)qΔijkxx=Nixqi′x+Njxqj′x+Nkxqk′x,qΔijkyx=Nixqi′y+Njxqj′y+Nkxqk′y,qΔijkzx=Nixqi′z+Njxqj′z+Nkxqk′z,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). The heat flux along the normal direction of the element face is calculated by(43)qΔijk=−qΔijkxnx+qΔijkyny+qΔijkznz,where nx, ny, and nz are the direction cosines of outward normal vector of the element face. The 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 Pi2e in equation (22) as(44)PΔijkic=∫ΔijkNiqΔijkdsSince 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)(45)∂TΔijkic∂t=PΔijkicCii=∫ΔijkNiqΔijkdsCii.The above formula is given for nodei of the element face. The 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(46)∂Tic∂t=∑j=1nΔ∫ΔjNiqΔjdsCii,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(47)∂Ti∂t=1CiiPi1+Pi2+Pi3+Pi4+Pic,where(48)Pic=∑j=1nΔ∫ΔjNiqΔjds.The integral in the above equation is calculated by the three-point symmetric numerical integration formulae proposed by Hammer et al. [28].The above coupling algorithm is simple in implementation. Though FEM and the meshless particle method employ approximations of first-order consistency in this work, the coupling algorithm is applicable when higher-order 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).


## 4.2. 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(49)Tit+Δt=Tit+∂Tit∂tΔt,for FEM, and with(50)Tit+Δt=Tit+∂Tcit∂tΔt,for the meshless particle method, where Δt is the time step. The current work uses a fixed time step for time integration.The forward Euler method is conditionally stable, so the time step used for calculation must be less than a critical time step. The critical time step is estimated for three-dimensional problems as(51)Δtc=βmincρlEmin2λ,cρlPmin2λ,where β is a scale factor, lEmin is the minimum element size characterized by minimum altitude for tetrahedron element, and lPmin 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. The 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.


## 5. Numerical Validation

A numerical example concerning heat transfer in a solid cube is used to examine the correctness and accuracy of the coupled method. The solid cube has initial temperatureTx,0=1K everywhere. It is exposed to the surroundings with temperature Ta=0K. The 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=1W/m2⋅K. The cube has a side length of 2a (a = 1 m); its center is at the coordinate origin. Its material properties are assumed to be ρ=1kg/m3, c=1J/kg⋅K, and λ=1W/m⋅K. There 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. The detailed expression for the analytical solution is given as(52)Tx,t=φx,tφy,tφz,t,where(53)φx,t=∑i=1∞sinnianiaμicosnixe−kni2t.In equation (53), k=λ/ρc, and ni denotes the positive roots of the following equation:(54)natanna=haλ,μi=121+sin2nia2nia.φ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 Figure4(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. The meshless region is discretized with 1,728 (12 × 12 × 12) uniform particles. The particle spacing is about 0.09 m. The FEM region is discretized with 61,070 elements, which have an approximate global size of 0.09 m. The total simulation time is set to 0.4 s.Figure 4
Calculation models for heat transfer in a solid cube: (a) coupled model; (b) FEM model; (c) meshless model.
(a)(b)(c)Firstly, we choose the coefficient for smoothing lengthαh as 1.5 and the critical radius rcrit for determining ghost particles as twice smoothing length (about 0.27 m) for a preliminary correctness test of the coupled method. In this circumstance, lPmin and lEmin are 0.136 m and 0.022 m, respectively. The critical time step Δtc 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. This 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. The maximum relative error is defined as(55)est=max1≤i≤ntTnxi,t−Taxi,tTaxi,t,and the global relative error is defined as(56)egt=∑i=1ntTnxi,t−Taxi,t2∑i=1ntTaxi,t2,where Tnxi,t and Taxi,t are the numerical and analytical results of temperature at node/particle i, respectively, and nt is the total number of particles and nodes used in numerical calculation. Table 1 shows the calculated es and eg at different time. As a reference, results of FEM and the meshless particle method are also provided in the table. The FEM model and meshless model are given in Figures 4(b) and 4(c), respectively. The 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. The 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 eg and es for the coupled method are only 0.0039 and 0.0578, respectively. This 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. This surface belongs to the FEM region. The 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. This 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.Figure 5
Temperature time histories of two observation points.Table 1
Relative error at different time points.
Methodes0.1eg0.1es0.2eg0.2es0.3eg0.3es0.4eg0.4Coupled method0.02040.00210.03890.00260.04900.00320.05780.0039FEM0.00870.00110.01010.00090.01080.00090.01110.0009The meshless particle method (αh = 1.5)0.00770.00150.00810.00100.00790.00070.00810.0009Figure 6
Relative error of temperature att = 0.4 s using different critical radius rcrit (rcrit = αr times smoothing length): (a) top surface, αr = 2; (b) plane of z = 0.045 m, αr = 2; (c) top surface, αr = 4; (d) plane of z = 0.045 m, αr = 4.
(a)(b)(c)(d)Then, parametric effects on the accuracy of the coupled method are investigated to identify optimal parametric values. Effect of critical radiusrcrit is firstly studied. rcrit 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, es decreases obviously, and eg also decreases on the whole. This means that the accuracy of the coupled method is improved with the increase of αr. The 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. This 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. This 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.Table 2
Relative error at different time points for different critical distances.
αres0.1eg0.1es0.2eg0.2es0.3eg0.3es0.4eg0.420.02040.00210.03890.00260.04900.00320.05780.003930.00860.00240.01170.00210.01340.00170.01610.001640.00860.00250.01170.00220.01340.00160.01410.001550.00860.00250.01170.00220.01340.00160.01410.0015Effect 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. The 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.Table 3
Relative error at different time points for different time steps.
∆tes0.1eg0.1es0.2eg0.2es0.3eg0.3es0.4eg0.45.0 × 10−40.00920.00240.01210.00220.01370.00180.01440.00162.5 × 10−40.00860.00250.01170.00220.01340.00160.01410.0015Table 4
Relative error at different time points for different smoothing lengths.
αhes0.1eg0.1es0.2eg0.2es0.3eg0.3es0.4eg0.41.00.00920.00110.01180.00090.01290.00090.01340.00091.250.00920.00160.01190.00150.01330.00130.01390.00131.50.00920.00240.01210.00220.01370.00180.01440.0016


## 6. Transient Heat Transfer Modeling for FSW of Al 6061 Plates

### 6.1. 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). The two aluminum alloy plates have the same dimensions which are 254 mm (length) × 102 mm (width) × 12.7 mm (thickness). The tool is made of H-13 tool steel. The main geometric dimensions of the tool are 6 mm of pin radiusrP, 12 mm of pin height hP, and 25 mm of shoulder radius rsh. The 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 Section2, a half model is established to simulate the FSW process by making use of the symmetry. Figure 7(a) shows the coupled model for the problem. The region of the workpiece close to the moving path of the tool is discretized by 96,658 particles (called real particles later for convenience). This region has a width of 30 mm which is beyond the shoulder’s radius. The rest region of the workpiece is discretized by 110,115 elements which have an approximate global size of 2 mm. Real particles except those on the coupling interface have a spacing of about 1 mm. Ghost particles and the real particles on the coupling interface have a spacing about 2 mm. Smoothing lengths of all particles are taken as 1.25 times particle spacings, i.e., 1.25 mm for most real particles and 2.50 mm for ghost particles and the real particles on the coupling interfaces. The critical radius rcrit for determining ghost particles is taken as 4 times smoothing length of real particles on the coupling interface, i.e., 10 mm. lPmin and lEmin for the current coupled model are 0.85 mm and 1.25 mm, respectively. The critical time step Δtc is estimated as 8 × 10−3 s with the scale factor β taken as 1. The 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 Δtc. The time step for the current calculation is finally set to 5 × 10−3 s. The total simulation time is 200 s. For comparison, a meshless model is also established, as shown in Figure 7(b). The meshless particle method uses the same time step for calculation. The 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.Figure 7
Calculation models: (a) coupled model; (b) meshless model.
(a)(b)Constants of heat source model used in the current calculation are given in Table5. The workpiece’s density is 2,700 kg/m3. Temperature dependence of specific heat, thermal conductivity, and yield stress are considered for the workpiece material. The specific heat and thermal conductivity are defined as a cubic function of temperature. The temperature dependence of yield stress is defined in a tabulated form. The parameters in temperature-dependent specific heat, thermal conductivity, and yield stress are taken with the same values as those in Xiao et al. [12]. The convective heat transfer coefficients h1 and h2 are 30 W/(m2·K) and 200 W/(m2·K), respectively.Table 5
Heat source model constants.
Parameter namets (mm)tb (mm)CfParameter value220.85


### 6.2. Results and Discussion

Figure8 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. Thus, 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. Thus, the top region of the workpiece within the shoulder becomes the highest temperature zone. As shown in Figures 8(c)–8(f), temperature distribution in the welding region gradually stabilizes with the increase of welding time. Temperature distributions shown in Figures 8(e) and 8(f) are quite similar in the welding region. This 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. The 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]. They also agree well with those presented in the current work.Figure 8
Temperature contours for the results of the coupled method: (a)t = 2 s; (b) t = 3 s; (c) t = 10 s; (d) t = 20 s; (e) t = 30 s; (f): t = 40 s.
(a)(b)(c)(d)(e)(f)Figure 9
Temperature contours for the results of the meshless particle method: (a)t = 2 s; (b) t = 3 s; (c) t = 10 s; (d) t = 20 s; (e) t = 30 s; (f) t = 40 s.
(a)(b)(c)(d)(e)(f)Figure10 shows thermal cycles obtained by calculation and experiment for four typical points. Points 1 and 2 are both 2 mm under the workpiece’s top surface, and they are 8 mm and 16 mm from the weld line, respectively. Points 3 and 4 are both 8 mm under the workpiece’s top surface, and they are 8 mm and 16 mm from the weld line, respectively. As seen from Figure 10, the results calculated by the coupled method are in good agreement with the experimental results. The 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 10
Thermal cycles for the four points: (a) Point 1; (b) Point 2; (c) Point 3; (d) Point 4.
(a)(b)(c)(d)Figure11 gives the change of maximum temperature of the workpiece. The result of the coupled method is in good agreement with the meshless particle method. The maximum temperature rises sharply in the tool plunge process and the initial stage of welding process. Then, 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. The 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.Figure 11
Change of maximum temperature of the workpiece.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. The 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. The 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. Thus, the coupled method will be more competitive over the meshless particle method for thermal-mechanical analysis.


## 6.1. 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). The two aluminum alloy plates have the same dimensions which are 254 mm (length) × 102 mm (width) × 12.7 mm (thickness). The tool is made of H-13 tool steel. The main geometric dimensions of the tool are 6 mm of pin radiusrP, 12 mm of pin height hP, and 25 mm of shoulder radius rsh. The 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 Section2, a half model is established to simulate the FSW process by making use of the symmetry. Figure 7(a) shows the coupled model for the problem. The region of the workpiece close to the moving path of the tool is discretized by 96,658 particles (called real particles later for convenience). This region has a width of 30 mm which is beyond the shoulder’s radius. The rest region of the workpiece is discretized by 110,115 elements which have an approximate global size of 2 mm. Real particles except those on the coupling interface have a spacing of about 1 mm. Ghost particles and the real particles on the coupling interface have a spacing about 2 mm. Smoothing lengths of all particles are taken as 1.25 times particle spacings, i.e., 1.25 mm for most real particles and 2.50 mm for ghost particles and the real particles on the coupling interfaces. The critical radius rcrit for determining ghost particles is taken as 4 times smoothing length of real particles on the coupling interface, i.e., 10 mm. lPmin and lEmin for the current coupled model are 0.85 mm and 1.25 mm, respectively. The critical time step Δtc is estimated as 8 × 10−3 s with the scale factor β taken as 1. The 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 Δtc. The time step for the current calculation is finally set to 5 × 10−3 s. The total simulation time is 200 s. For comparison, a meshless model is also established, as shown in Figure 7(b). The meshless particle method uses the same time step for calculation. The 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.Figure 7
Calculation models: (a) coupled model; (b) meshless model.
(a)(b)Constants of heat source model used in the current calculation are given in Table5. The workpiece’s density is 2,700 kg/m3. Temperature dependence of specific heat, thermal conductivity, and yield stress are considered for the workpiece material. The specific heat and thermal conductivity are defined as a cubic function of temperature. The temperature dependence of yield stress is defined in a tabulated form. The parameters in temperature-dependent specific heat, thermal conductivity, and yield stress are taken with the same values as those in Xiao et al. [12]. The convective heat transfer coefficients h1 and h2 are 30 W/(m2·K) and 200 W/(m2·K), respectively.Table 5
Heat source model constants.
Parameter namets (mm)tb (mm)CfParameter value220.85


## 6.2. Results and Discussion

Figure8 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. Thus, 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. Thus, the top region of the workpiece within the shoulder becomes the highest temperature zone. As shown in Figures 8(c)–8(f), temperature distribution in the welding region gradually stabilizes with the increase of welding time. Temperature distributions shown in Figures 8(e) and 8(f) are quite similar in the welding region. This 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. The 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]. They also agree well with those presented in the current work.Figure 8
Temperature contours for the results of the coupled method: (a)t = 2 s; (b) t = 3 s; (c) t = 10 s; (d) t = 20 s; (e) t = 30 s; (f): t = 40 s.
(a)(b)(c)(d)(e)(f)Figure 9
Temperature contours for the results of the meshless particle method: (a)t = 2 s; (b) t = 3 s; (c) t = 10 s; (d) t = 20 s; (e) t = 30 s; (f) t = 40 s.
(a)(b)(c)(d)(e)(f)Figure10 shows thermal cycles obtained by calculation and experiment for four typical points. Points 1 and 2 are both 2 mm under the workpiece’s top surface, and they are 8 mm and 16 mm from the weld line, respectively. Points 3 and 4 are both 8 mm under the workpiece’s top surface, and they are 8 mm and 16 mm from the weld line, respectively. As seen from Figure 10, the results calculated by the coupled method are in good agreement with the experimental results. The 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 10
Thermal cycles for the four points: (a) Point 1; (b) Point 2; (c) Point 3; (d) Point 4.
(a)(b)(c)(d)Figure11 gives the change of maximum temperature of the workpiece. The result of the coupled method is in good agreement with the meshless particle method. The maximum temperature rises sharply in the tool plunge process and the initial stage of welding process. Then, 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. The 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.Figure 11
Change of maximum temperature of the workpiece.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. The 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. The 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. Thus, the coupled method will be more competitive over the meshless particle method for thermal-mechanical analysis.


## 7. Conclusions

A coupled method of FEM and meshless particle method is developed for simulating heat transfer process during FSW. The 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. The study of parametric effect shows that its accuracy is improved with the increase of critical radiusrcrit for determining ghost particles till it reaches about 4 times smoothing length, and has no significant dependence on time step. The 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.The coupled method is successfully applied to the simulation of heat transfer during FSW of Al 6061-T6 plates. Thermal 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.The 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.



---

*Source: 2574127-2020-05-20.xml*