An alternative pressure-dependent velocity boundary condition for modeling self-reacting friction stir welding

A pressure-dependent velocity boundary condition is developed based on wear theory for modeling self-reacting friction stir welding using computational fluid dynamics approach, which provides a new perspective in understanding the physics of sliding/sticking transition condition. The importance of shear layer in weld formation is emphasized. Effects of welding speed on the weld cross-section geometry can be robustly captured with this newly developed boundary condition. Computational results show that at higher welding speed, the TMAZ boundary moves towards the pin periphery at the advancing side, which corresponds to the experimental observations. This tendency could serve as a numerical criterion to predict void defect formation.


Introduction
Self-reacting friction stir welding (SRFSW) is an advanced variant of friction stir welding (FSW), which utilizes a double-sided tool and shows several superiorities. The additional root shoulder can effectively avoid root penetration defects in conventional FSW welds. Meanwhile, the axial force between the crown and root shoulder surfaces is balanced, which results in a trivial net axial force and enables in-space FSW without requirement of heavy manufacturing equipment. SRFSW can be divided into two categories based on tool designs. First is the fixed gap bobbin tool [1], where the distance between the crown and root shoulders is constant. The second is the adaptive tool [2], where the distance between the crown and root shoulders can be dynamically controlled according to the force acting between the two shoulders. One typical application of SRFSW is manufacturing large-scale cryogenic fuel tanks in space vehicles.
Recognizing the importance of understanding the weld formation mechanism during FSW/SRFSW, various researches have been taken to investigate the experimentally observed shear layer [3][4][5] formed in the process. As a result of the interactions between the tool and the workpiece, localized shear develops in the tool vicinity and forms a narrow region of intense plastic deformation, i.e., shear layer, which is also referred to as shear band [6,7]. Shear layer formation in sheardominated processes is generally believed to be first observed by Tresca [8]. To understand this deformation localization mechanism, Batra and Wei [6] analyzed the initiation and development of shear layer using one-dimensional (1D) transient thermo-elasto-viscoplastic formulation, with prescribed tangential velocity and constant heat flux as boundary conditions. They found that the shear band initiation time closely depends upon the prescribed heat flux. Pei and Dong [9] incorporated this 1D shear localization model into FSW process modeling, with different boundary and initial conditions. In their following study [10,11], they used Zener-Hollomon material constitutive model instead of the Johnson-Cook model, which was previously used in [6,9]. They calculated the width and formation time of shear layer, which can be served as criterions for quality weld.
Compared with the shear localization model which only provides information surrounding the tool, thermomechanical computational fluid dynamics (CFD) models are capable of predicting the entire weld cross-section morphology [12,13], visualizing flow field [13][14][15][16][17][18], and obtaining weld thermal history [13-15, 19, 20]. Due to the complexities of the sticking/sliding transition behavior between tool and workpiece, different assumptions and simplifications have been made at tool/workpiece (T/W) boundary. One is velocity boundary condition, where a prescribed velocity is assigned at T/W interface. The value can be equal to the tool velocity [14,21] representing a full sticking condition, or a fraction of tool velocity representing a partially sticking condition. The fraction is either a constant [22] or a preassigned variable as an (empirical) function of the radial distance [23][24][25][26]. The other group of T/W boundary condition is based on friction stress [18,20,27], and the velocity of interfacial material is solved implicitly. Chen et al. [15] proposed a relative velocityconstrained shear stress boundary condition. Andrade et al. [18] investigated the effect of parameters of a parametrical relative velocity-constrained shear stress boundary condition by evaluating the velocity distribution at the T/W interface. Zhao and Liu [13] evaluated the credibility of velocity and shear stress boundary condition, respectively, based on the calculated weld cross-sectional morphology. According to their research, the shear stress boundary condition represents T/W sliding in Coulomb friction and subsurface interfacial sliding in Tresca friction, whereas the velocity boundary condition represents an averaging estimated sticking/sliding transition state of the shear layer, which is more appropriate to model the actual T/W contact behavior. The importance of shear layer in quality weld formation is emphasized in their work. In addition to CFD models, computational solid mechanics (CSM) based thermo-mechanical models are of interest for researchers [28][29][30][31]. Dialami et al. [32] developed a two-stage solution strategy to solve the coupled thermomechanical problem. However, owing to the lack of fundamental understanding of the T/W friction behavior and challenges of direct experimental observations, an explainable velocity boundary condition which combines pressure distribution that can be implemented under various welding parameters and produce robust predictions in CFD-based SRFSW/ FSW process models is still missing.
The focus of this study is on developing a new pressuredependent velocity boundary condition based on wear theory, which provides a new and explainable perspective in understanding the physics of sliding/sticking condition and shear layer formation mechanism during FSW/SRFSW CFD modeling. The corresponding weld cross-sectional morphology is predicted based on plastic strain distribution and validated with experimental macrographs at different welding speeds. The rationale behind pure velocity and pressure-dependent velocity boundary condition will be compared and evaluated.

Governing equations
In CFD-based FSW/SRFSW modeling, solid-state metals are treated as non-Newtonian fluids with high viscosities. For incompressible single-phase fluid, governing equations of the mass, momentum, and energy conservation equations are as follows: where v ! is the velocity vector, ρ the density, t the time, p the pressure, C p the specific heat capacity, k the thermal conductivity, and S v the spatial heat source term from plastic deformation. μ is the effective viscosity. The viscoplastic material mechanical behavior is described by the Sellars-McTegart constitutive model [33]: where σ e is the flow stress, Q the deformation activation energy, R the gas constant, and T the absolute temperature; A, α, and n are material constants. Effective strain rate is given by [21,34]: where D is the strain rate tensor: In incompressible single-phase flow, ∂v i ∂x i ¼ 0. Effective viscosity μ could be calculated from effective stress σ e and effective strain rateε e as [21,27,35] Flow stress in Eq. (4) could also be rewritten as [36] where Z is the Zener-Hollomon parameter, which represents the temperature compensated effective strain rate and is expressed as Material parameters in constitutive models and thermophysical properties of AA6061-T6 used in this research are the same as in [13] and are according to tensile tests data in [20]. Constitutive parameters of AA2219-T87 implemented in this paper are the same as in [37]. Density, heat capacity, and thermal conductivity of AA2219-T87 are taken as 2840kg· m −3 , 864J·Kg −1 ·K −1 , and 121W·m −1 ·K −1 , respectively. An artificial term is added in Eq. (8) to enforce flow stress decreases to zero as temperature approaches the solidus.

Boundary conditions
In this study, two mechanical T/W boundary conditions are investigated. The first one is the constant pure velocity boundary condition, and the second is the pressure-dependent velocity boundary condition, which is derived from the modified Archard's wear equation.
In the constant velocity boundary, the material velocity at interface is directly assigned as a fraction of tool velocity as where δ represents the degree of sticking and an averaging effect of the shear layer [13]. As δ increases, the predicted width of the deformation zone increases accordingly since higher boundary velocity intensifies the material movement in the pin vicinity. By adjusting modeling results to match with the experimental weld cross-section morphology in [38], the value of δ is determined as 0.013. The second T/W boundary condition is based on wear theory. Conventionally, wear models are applied in FSW to predict FSW tool wear [39]. In contrary, this study adopts an opposite perspective and considers the FSW/SRFSW process as the wear of the workpiece. In other words, to capture the effects of shear layer, the wear model is implemented on the workpiece instead of the tool. This assumption is reasonable since the tool always has higher hardness than the workpiece material. As the workpiece materials are worn by the friction action of the tool, the worn metal forms a rotating shear layer surrounding the tool and deposited backward in the weld. Figure 2 shows the experimental observations of shear layer during FSW. The importance of this shear layer formation to achieve defect-free welds has also been discussed in [13].
Wear damage of the workpiece can be described by Archard's equation [40]: where V is the total wear volume; k is the dimensional Archard wear coefficient, which depends on the material and surface cleanliness; W is the applied load; and S is the sliding distance between the tool and the workpiece. The thickness of shear layer can be represented by the wear depth. Therefore, the modified Archard's equation [39] is adopted as: where h is the wear depth, p the contact pressure, t the time, and v slip the relative velocity between the tool and the workpiece. According to Eq. (11) or Eq. (12), in Archard's wear model, the worn material/shear layer is continuously generated and has the same velocity as the tool, which represents the full sticking scenario. However, during FSW/SRFSW, it is known that the T/W interface is subject to a complex sticking/sliding transition state, where the sticking and sliding condition varies during each tool revolution at tool surface [13]. As shown in Figure 1, in the scenario of full sticking, the velocity of the material at T/W interface equals the tool velocity, and the shear layer continuously forms. When pure sliding happens, the wear of workpiece does not occur, and the generation of the shear layer stops, and the velocity of the formed shear layer equals to zero. From the modeling perspective, the velocity of the shear layer is supposed to be estimated and applied as the T/W workpiece boundary condition. Under the same rotation speed of the tool, the sliding distance between the shear layer and the remaining material of the workpiece can be written as and the relative velocity between the shear layer and the remaining material can be obtained as where β(p, t) defined as wear fraction coefficient, which is a function of pressure and time, represents the fraction of wear that happened along the tool traveled distance. βv tool represents the velocity of the shear layer v shear layer . S tool is the Fig. 1 Illustration of sticking/sliding transition state at T/W interface: a full sticking; b pure sliding distance the tool moves; S mat is the distance remaining material moves; v tool is the tool velocity; v mat is the remaining material velocity. Different from a constant value of δ in Eq. (10), β(p, t) is dependent on the wear condition, which should be a function of pressure and time based on Archard's wear model. β(p, t) is an advanced version of δ and provides a more accurate estimation of the equivalent shear layer velocity in sticking/sliding transition state. β(p, t) and δ have the same magnitude.
Since the relative velocity is time-dependent, furthermore, the SRFSW is treated as a quasi-steady state in process modeling, the wear depth can be calculated as where t 1 is the moment when shear layer starts to form and t 2 is the moment when the shear layer detaches from the tool. β p ð Þ is the time-averaged wear fraction coefficient of β(p, t), which is a function of pressure; v mat is the averaged remaining material velocity; Δt = t 2 − t 1 is the total existing time of each shear layer. v mat and Δt are constants under the same welding condition.
During FSW, a newly formed shear layer is established and circulates around the pin for a short distance to the trailing side until it deposits [4]. This is defined as one wear cycle, which happens every tool revolution. When the welding speed is constant, the average thickness of the worn material during each wear cycle remains constant to form a weld. This agrees well with experimental observations of the uniform thickness of shear layer during each wear cycle [4,5,41] as shown in Figure 2.
When the welding speed varies within a reasonable range, based on experimental observations from literature [38,42], void defects always occur at the advancing side due to the lack of deposition, i.e., the shear layer could not be transported to a sufficient distance into the advancing side, rather than occur between shear layers. This indicates the amount of shear layer generated is adequate under different welding speeds without leaving lack of fills in between successive layers. Therefore, according to the geometrical restriction, at the same rotation speed, the thickness of the shear layer could be assumed to be proportional to the welding speed [4,43].
Therefore, according to Eq. (15), the pressure-dependent velocity boundary condition can be formulated as where α and γ are constants and v welding is the welding speed. Equation (16) only holds when p is larger than a threshold value. Formation of shear layer would not occur when p is low, and the boundary velocity would decrease to zero accordingly. The values of α and γ are determined by calibrating the modeling results with the experimental weld cross-section morphology. It should be noted that continuity is required in the CFD computational domain. Accordingly, both pure velocity boundary condition and pressure-dependent velocity boundary condition represent an estimation of the equivalent shear layer velocity to capture the effects of the shear layer on the remaining material. The shear layer itself is not explicitly modeled. From thermal perspective, frictional heat generation rate at T/W interface is calculated as where ζ is the heat partition coefficient, which represents the fraction of heat absorbed by the workpiece compared with total amount of generated heat. τ b is the shear stress based on Coulumb-Tresca friction model as decribed in reference [13]. ζ can be calculated by [24] ζ ¼ Fig. 2 a Top view of shear layer in FSW [41]. b shear layer behind the pin in FSW [5]. c Schematic illustration of shear layer around the pin [4] The heat partition coefficient at the bottom shoulder is 0.9. A slightly smaller heat partition coefficient of 0.8 is selected at the top shoulder since it is clamped into the entire machine, which acts as a large heat sink leading to higher amount of heat dissipation.
On the remaining surfaces of workpiece, convective heat transfer boundary conditions are adopted. The convection heat transfer coefficient and environment temperature are set as 20W/m 2 K and 300K, respectively.

Numerical model geometry and particle tracing methodology
A schematic illustration of the computational domain for the AA6061-T6 plate is provided in Figure 3. The SRFSW tool is hollowed out from the workpiece. Dimensions and geometry of the tool and the workpiece are consistent with the experimental configurations described in [38]. The calculation domain contains an inlet and an outlet surface indicating the direction of welding speed. Thickness of the workpiece is 4mm. The diameters of the non-threaded pin and top and bottom shoulder are 8mm, 18mm, and 16mm, respectively. Plunge depth of both top and bottom shoulders is 0.1mm. To validate the feasibility of the pressure-dependent velocity boundary, another set of experimental data for AA2219-T87 welds with a thickness of 6mm [44] is verified. The structure of the model geometry is similar to Figure 3. The diameter of the pin is 8mm, and both the top and bottom shoulders have a diameter of 16mm. Plunge depth of both top and bottom shoulders is 0.1mm. There are three flats on the smooth cylindrical pin. In this research, a cylindrical pin tool is hollowed out from the workpiece to represent the dynamic volume. To capture the intensified material velocity around the pin caused by the cutting effect from the flats [45], α in Eq. (16) is magnified, which will be discussed in the following section.
Interactions between tool and workpiece, which are incorporated with corresponding thermo-mechanical boundary conditions and heat generation rate described in Eq. (10), Eq. (16), and Eq. (17), as well as material properties described in Eq. (7), are encoded with user-defined functions (UDF) and implemented into the FLUENT program. In this study, four cases for the AA6061-T6 plate and two cases for the AA2219-T87 plate are calculated and compared, as summarized in Table 1.
Particle tracing methodology is adopted according to reference [13] to determine the strain distribution on the weld section behind the pin, after the simulation entered the quasi-steady state. A virtual receiving plane of 20mm × 3.8mm is placed at 30mm behind the tool center, where the stirring effects of pin on material flow diminished. A total number of 1600 particles are traced. The strain can be calculated by integrating the effective strain rate over time. By combining the relationship between time and velocity, the strain at each point on the receiving plane can be integrated along the path lines as where v is the local velocity magnitude.

Results and discussion
3.1 Strain distribution and weld cross-sectional morphology prediction Figures 4,5,6,7,8,and 9 show the strain distribution on the weld cross-section at 30mm behind the tool center of the six simulated cases. The tool profile is marked out with red dotted lines on the calculated iso-strain contour maps and experimental weld macrographs. Boundary of deformation zone, including both stirring zone (SZ) and thermo-mechanically affected zone (TMAZ), can be sketched with iso-strain contours since these two zones were subject to mechanical deformation during the welding process. At the advancing side, the iso-strain Fig. 3 Illustration of model geometry contours are denser than the retreating side, indicating higher spatial gradient of strain. This is consistent with the general SRFSW weld macrograph observations, where the boundary between TMAZ and HAZ is more definable in the advancing side whereas gradual transition occurs at the retreating side. Comparing Figures 4 and 5, under the same welding speed of 100mm/min, pure velocity boundary condition (case 1) and pressure-dependent velocity boundary condition (case 2) predict similar iso-strain contours at the advancing side, both of which are in line with the experimentally observed TMAZ boundary. The distance between the pin periphery and the iso-strain contour at the neck of the hourglass morphology is also in agreement with the weld macrograph. In contrary, the morphology in the retreating side is hardly captured with the pure velocity boundary, where the high strain contours deviate significantly in the crown shoulder region, as shown on the left side of Figure 4. This is greatly improved with the newly developed pressure-dependent boundary condition, as shown in Figures 5 and 8.
As welding speed increases from 100mm/min to 250mm/min, based on the pressure-dependent velocity boundary, the predicted TMAZ boundary in the advancing side moves towards the pin center, as shown from Figures 5, 6, and 7. This trend is consistent with the experimental observations. At the welding speed of 250mm/min in case 4, the TMAZ boundary lies within the pin periphery in both the computational results and the experimental weld macrographs, as shown in Figure 7. This high welding speed parameter actually leads to the formation of void defect, which is located within the pin periphery in the advancing side as in Figure 7c. Generally, in the SRFSW process, the welding parameters that lead to void defects formed at this location would also result in the TMAZ boundary lies a close distance within the pin periphery in the advancing side, as shown in Figure 10. Accordingly, the relative position of the predicted TMAZ boundary to the pin periphery in the advancing side could serve as an effective criterion for void defect prediction.  Fig. 4 a Cloud map of predicted strain distribution at 30mm behind tool center in case 1. b Isostrain contour map at 30mm behind tool center in case 1. c Crosssection of the weld under welding speed of 100mm/min from reference [38] For AA2219-T87 welds, as shown in Figure 9b and d, the predicted TMAZ boundary at advancing side shifts closer to the pin center when welding speed increases from 150 to 350mm/min, which agrees with the experimental weld cross-section morphology in Figure 9a and c. This indicates good applicability of the pressure-dependent velocity boundary when different material and welding parameters are adopted.
As discussed in Section 2.2, the values of α and γ are determined by calibrating the modeling results with the experimental weld cross-section morphology. In case 1 to case 4, the values of α and γ are 9.5 × 10 −3 and 4.5 × 10 5 when the unit of v welding is mm/min. In case 5 and case 6, they are 2.684 × 10 −2 and 2.3 × 10 5 , respectively. Since the material velocity around the pin is intensified by the flats, compared with the value in case 1 to case 4 where a smooth cylindrical Fig. 5 a Cloud map of predicted strain distribution at 30mm behind tool center in case 2. b Isostrain contour map at 30mm behind tool center in case 2. c Crosssection of the weld under welding speed of 100mm/min from reference [38] Fig. 6 a Cloud map of predicted strain distribution at 30mm behind tool center in case 3. b Isostrain contour map at 30mm behind tool center in case 3. c Crosssection of the weld under welding speed of 150mm/min from reference [38] pin is adopted, α is increased by more than 2.5 times to capture the effects of the flats. This indicates in the CFD model, the equivalent effects of the flats on the material flow field can be included as an intensified velocity boundary condition. Fig. 7 a Cloud map of predicted strain distribution at 30mm behind tool center in case 4. b Isostrain contour map at 30mm behind tool center in case 4. c Crosssection of the weld under welding speed of 250mm/min from reference [38] Fig. 8 a Cloud map of predicted strain distribution at 30mm behind tool center in case 5. b Isostrain contour map at 30mm behind tool center in case 5. c Crosssection of the weld under welding speed of 150mm/min from reference [44]

Flow field and pressure distribution
The material flow fields for quality welds are overall similar. To evaluate the effectiveness of the pressuredependent boundary condition in capturing the results of different welding speeds, especially the difference when void defect forms, case 1 to case 4 will be focused on and discussed. Material flow fields of case 1 to case 4 are provided in Figure 11. The path line distributions are overall similar. Incoming materials are mainly plasticized at the retreating side and some fractions are deposited backward at the advancing side. The red arrow indicates the position where the materials start to deposit. As the welding speed increases, the deposition position moves away from the shoulder edge in the advancing side towards the pin center, which means less amount of material is transported to a sufficient distance to fill up the advancing side. This tendency matches well with the weld cross-section morphology from Figures 5, 6, and 7. Since less amount of material is stirred into the advancing side, the deformation zone size is narrowed and leads to a TMAZ boundary lies within the pin periphery. When the amount of material filled into advancing side drops below a certain degree, void defect is formed as shown in Figures 7 and 11d.
The velocity distribution on the mid-pin depth plane is shown in Figure 12. As the welding speed increases from Figure 12b to d, overall, the magnitude of material velocity increases. The red arrow in these figures points towards the low-velocity region and its rotating direction agrees with that of flow field distribution (red arrow in Figure 11) as the welding speed increases. Under pure velocity boundary condition in Figure 12a, the magnitude of velocity is distributed Fig. 9 a Cross-section of the weld under welding speed of 150mm/min from reference [44]. b Iso-strain contour map of predicted strain distribution at 30mm behind tool center in case 5. c Cross-section of the weld under welding speed of 350mm/min from reference [44]. d Iso-strain contour map of predicted strain distribution at 30mm behind tool center in case 6 Fig. 10 Cross-section of the welds of AA6061-T6 from reference [42] relatively symmetric in front and behind the pin. With the pressure-dependent velocity boundary condition, the shape of velocity magnitude field is distorted and points towards the position where materials start to deposit at the advancing side, as discussed in Figure 11. Comparing Figure 12a and b, under the same welding speed of 100mm/min, the maximum  Figure 13 shows the temperature distribution in the center of the tool cross-section using pressure-dependent velocity boundary condition at the welding speed of 100mm/min. Thermal histories in the middle thickness of the workpiece are extracted at two points located at pin peripheries on both advancing and retreating sides, as illustrated in Figure 13. Figure 14 compares the corresponding thermal histories using pure velocity boundary and pressure-dependent velocity boundary under the welding speed of 100mm/min. At the pin periphery, the maximum temperature approaches the solidus of AA6061-T6 at either the advancing or retreating side. With the velocity boundary condition, temperature in the advancing side is higher than that in the retreating side around the pin periphery. In contrary, with the pressure-dependent velocity boundary condition, the retreating side has a higher temperature. The temperature results are closely related to the material flow behavior. With the pressure-dependent velocity boundary condition, the existence of low-pressure region limits the amount of material flow towards the advancing side, as shown in the flow field in Figure 11b. Based on the modeling results of the flow field, plasticized material is mainly deformed at the retreating side and result in a higher temperature on the retreating side. With the pure velocity boundary condition, material flow is enforced in the vicinity of the pin and is transported further into the advancing side. Materials experience longer period of friction and deformation heating, which results in higher temperature at the advancing side.

Conclusions
In this paper, a new pressure-dependent velocity boundary condition based on Archard's wear equation is developed to model SRFSW using the CFD method. Derivation of this new T/W boundary condition is from the perspective of shear layer formation and the sticking/sliding transition state at the T/W interface is defined through the wear fraction coefficient. The main conclusions are as follows: & The weld cross-section geometry can be successfully predicted based on the strain distribution calculated from the pressure-dependent velocity boundary condition. The hourglass shape and location of the TMAZ boundary at the advancing side match well with the experimental results. & The pressure-dependent velocity boundary condition can robustly capture the effects of different welding speeds ranging from 100 to 250mm/min on the weld geometry. The relative position of the calculated TMAZ boundary to the pin periphery at the advancing side could serve as a numerical criterion to predict void defect formation. & Compared with pure velocity boundary condition, the pressure-dependent velocity boundary condition explicitly considers the varying pressure surrounding the pin and leads to more reasonable computed results of pressure and flow field distribution. & In CFD models, the equivalent effects of the flats on the material flow field can be included as an intensified velocity boundary condition.
More experiments will be performed to further validate the modeling results under different welding parameters and provide accurate quantitative predictions.