An integrated multiphysics model for friction stir welding of 6061 Aluminum alloy

This article presents a new, combined ‘integrated’‘multiphysics’ model of friction stir welding (FSW) where a set of governing equations from nonNewtonian incompressible fluid dynamics, conductive and convective heat transfer, and plain stress solid mechanics have been coupled for calculating the process variables and material behaviour both during and after welding. More specifically, regarding the multiphysics feature, the model is capable of simultaneously predicting the local distribution, location and magnitude of maximum temperature, strain, and strain rate fields around the tool pin during the process; while for the integrated (post-analysis) part, the above predictions have been used to study the microstructure and residual stress field of welded parts within the same developed code. A slip/stick condition between the tool and workpiece, friction and deformation heat source, convection and conduction heat transfer in the workpiece, a solid mechanics-based viscosity definition, and the Zener-Hollomonbased rigid-viscoplastic material properties with solidus cut-off temperature and empirical softening regime have been employed. In order to validate all the predicted variables collectively, the model has been compared to a series of published case studies on individual/limited set of variables, as well as in-house experiments on FSW of aluminum 6061.


INTRODUCTION
Friction Stir Welding (FSW), a solid state joining method developed and patented by TWI Ltd., Cambridge, UK in 1991 [1], has attracted significant interest from aircraft and car manufacturers for joining high strength aluminum alloy components.Specific examples of target materials for FSW include the wrought 6000-series Al-Mg-Si (Cu) alloys that are commonly used in aircraft fuselage skin and automotive inner body panels, mainly due to their ability to be strengthened by artificial aging after forming.FSW has also been used to produce rocket shells, the panel of the cabin of aircrafts with stringers and beams, hollow panels of wagons, and pipes [2].
Various numerical methods including the finite element, finite volume, finite difference, and analytical methods are utilized in developing prediction models of FSW.The governing equations (analytical or coupled with process parameters) are based on thermal [3], fluid dynamics [4][5][6][7], thermomechnical or solid mechanics [8][9][10][11][12][13][14][15], conceptual [16], and kinematics [17] models, often developed in-house or built in conventional codes such as Abaqus, COMSOL, Fluent, Deform-3d, CTH, Forge3, Ansys, Nisa, Star-ccm, iStir, Stir3d, Fidap, Sysweld, Cosmos, Hickory, Thercast, Fastflo, LS-Dyna, and Weldsim.The solid mechanics models consider that the welding material is solid during FSW and strain distribution is one of the model results.The fluid dynamics models normally consider a solid mechanics-based definition of viscosity and non-Newtonian incompressible flow of the softened solid.Hence, they can predict strain rate distributions and not the strain distribution; however there have also been some more other models which have used post processing techniques to compute strain on streamlines.Namely, Long et al [6] used a simple geometry-based formulation to calculate engineering strains on different streamlines and Arora et al. [7] estimated the accumulated strains experienced in the weld material by integrating strain rates over time along limited streamlines.In the multiphysics model in this article, based on a solid mechanics definition of viscosity we applied strain rate integration over time to compute plastic strains in all points of fluid dynamics model for the first time (not limited to few streamlines).
The assumptions made during the simulation of FSW can affect the analysis results considerably.Some of the common assumptions are: 2D or 3D modeling, assigned values of the process parameters, coupled/semi-coupled or thermal-pseudo-mechanical relations between thermal and mechanical properties, Lagrangian/Eulerian or Arbitrary Lagrangian Eulerian (ALE) formulations, transient or steady state conditions, consideration of the welding material as a fluid or a solid, assigning the governing equations by means of rigid viscoplastic or elastic-viscoplastic material models (such as Zener-Hollomon, Johnson-Cook, Norton-Hoff, etc.) for the solid state or viscosity models (such as power law, Carreau, Perzyna, etc.) for the fluid state, temperature-dependent thermal properties of the material, tool/workpiece interface conditions (slip/stick) [18], and workpiece boundary conditions.
There exist other modeling approaches such as 'integrated modeling' [19][20][21][22][23][24][25][26] and 'multiphysics modeling' [27][28][29][30][31][32][33] that have been applied in more recent investigations for modeling forming processes.In the latter, it has become more and more evident that the use of computer simulation is necessary for arriving at optimized products and structures with a minimum cost and time.Today, complex thermal manufacturing processes such as casting and welding are often studied with multiphysics models involving both computational fluid dynamics (CFD) and computational solid mechanics (CSM), as well as thermodynamic and kinetic models.In studies on welding processes, multiphysics modelling has also been combined with optimization methods to achieve desired properties of the final welds.In addition, the integration of multi-process steps and materials modeling has emerged as a new growing field [34][35][36][37][38][39].
The authors recently reported some advantages of the integrated multiphysics modeling in directly predicting microstructure and residual stress after FSW based on the main process parameters (weld speed and rotational speed) [40,41].There are two different regimes of material flow during FSW; namely "pin-driven flow" and "shoulder-driven flow".These material flow regimes should properly merge together to form a defect-free weld.In order to facilitate understanding the effect of parameters controlling these regimes during FSW, the main goal of the current article is to present other results of the 'integrated' 'multiphysics' model for simultaneously studying a multitude of variables of interest for a typical aluminum alloy.In multiphysics modeling, one uses different governing equations to predict the response variables during the forming process (such as temperature, strain, shear strain rate, etc); whereas in an integrated model, one uses the above main parameters to predict other (post-welding) parameters such as the weld microstructure and residual stresses [25].The present model brings these two prediction routines together under a single comprehensive numerical analysis framework, thus facilitating the study of any combination of main and post-process parameters.To show the applicability of the model, in addition to comparing to in-house experimental measurements, prediction results are tested against a set of case studies from the literature focusing on temperature, strain, and strain rate fields.

MODEL DESCRIPTION 2.1 GENERAL
The basic concept behind FSW is simple: A non-consumable rotating tool with a specially designed pin and shoulder is inserted into the abutting edges of the two parts to be joined and traversed along the line of joint (Fig. 1).When the plate is thick enough we may approximate a two dimensional flow in mid thickness of the plate and consider the shoulder as a remote heat source.The steady-state solution is a snapshot of the flow at an instant.

CFD GOVERNING EQUATIONS
To define the fluid-like behavior of the welding material around the pin, a set of governing equations assuming an incompressible flow has been used in the model as follows [42]: Momentum equation (Navier-Stokes equation): Dynamic viscosity: where u is the velocity vector, ρ is the density, p is pressure, F is the volume force vector, τ - is the effective shear stress, is the effective shear strain rate, and γ • is the dynamic viscosity.
The dynamic viscosity in Eq. (3), as opposed to a Newtonian fluid, may not be constant and can change by the shear flow stress and/or the shear strain rate (i.e., a non-Newtonian domain).This makes the coupled solution of partial deferential equations highly non-linear.As shown in Fig. 1, the pin is considered to rotate in its position (i.e., with no translational movement) and instead the plate moves from left to right.The inlet pressure at the left hand side of the plate and the outlet pressure at the right hand side of the plate are equal to ambient pressure as shown in Table 2.The upper and lower boundaries of the plate have moving-wall conditions equal to the weld transverse speed (u weld ).The outer boundaries of the plates are considered to be fixed and the material in contact with the pin is considered to accompany the pin during one rotation of the pin and then replaced by new material as the tool moves.This was based on a previous strain modeling study using experiments results via visioplastic strain measurements during FSW [43].It also complies with other evidence reported by Xu et al [44] that, in all cases involving a banded zone, the band thickness or spacing is found to be equal to the tool advancement during one pin rotation.
At the circularly shaped pin/plate interface (Fig. 2), a rotating velocity boundary condition is applied with the stick coefficient (δ) equal to 0.65 [45].The value was fitted by considering the torque measurement of the tool during FSW and its prediction by the model with different stick coefficients.The stick coefficient is initially defined as the ratio of tangential velocity of the material in contact with the pin to the tangential velocity of the pin periphery [5]: Considering a 2D model, the x and y components of the tangential velocity of the plate material are calculated using Eqs.( 5) and ( 6): Velocity in x direction: (5) Velocity in y direction: (6) where x and y are different positions of the pin based on the center of the coordinate system on the pin center, and ω is the pin rotational speed.

HEAT TRANSFER GOVERNING EQUATIONS
The governing differential equation of conductive and convective heat transfer in a fluid during steady state conditions is as follows.
Energy equation [46]: (7) where C p is the specific heat capacity at constant pressure, T is temperature, k is thermal conductivity, q is internal heat source.
According to the thermal boundary conditions shown in Fig. 1, the temperature at the left hand side of the plate is equal to ambient temperature (here 27 o C) and the right hand side of the plate has a convective heat flux temperature gradient.The upper and lower boundaries of the plate have insulation heat boundary conditions because of large dimensions of the plate compared to its thickness.The conductive heat transfer between the pin and plate is considered.The pin material is steel with a density of 7850 Kg/m 3 , specific heat capacity of 475 J/KgK, and thermal conductivity of 44.5 W/mK.Note that if heat transfer in the tool is not considered in the model, a lower value of temperature in the heading edge of the tool will be resulted.As in most pure thermal models of FSW [18], the heat generation from both plastic dissipation and frictional heat is modeled via a surface heat flux boundary condition at the tool-matrix interface via Eq.( 8).
Total heat flux [18]: (8) where r pin is the pin radius, ω is the tool rotational speed, δ is the stick coefficient, µ is the friction coefficient, and σ -is the effective normal stress (flow stress).

SOLID MECHANICS GOVERNING EQUATIONS
Governing equations of a rigid viscoplastic material at elevated temperature were used in the model via Eqs (9) to (13).
Zener-Hollomon equation [47,48]: Flow stress [49]: (10) Equivalent strain rate: (11) Strain rate tensor: (12) . ( ) ∇ = ∇ ∇ + Equivalent Shear stress: (13) where Z is the Zener-Hollomon parameter, ε • is the effective strain rate, Q is the temperature independent activation energy similar to that for self-diffusion, R is the gas constant, and α, A, and n are model constants determined by curve-fitting to hot deformation experimental data.Tello et al [50], for example, measured the flow stress in hot compression tests and determined the model constants of their materials.Here we also used an integration of strain rate over time to calculate plastic strain in all points of the model.
As stated earlier, friction stir welding is known to be a solid state welding process.This has been verified by earlier microstructural studies after the process where no evidence of dendrite microstructure was seen [44].Also there is no sharp decrease of transverse load during the process and hence the maximum temperature recorded during the process is normally below the solidus temperature of the material being welded.If we introduce a high volume of partial melting during FSW intentionally, it will generate a weld with low mechanical properties which is not desired.In order to implement the maximum temperature of solidus in the model, we consider that at a temperature equal or higher than solidus temperature, the flow stress (Eq.( 10)) is equal to zero, which in turn would cause no volumetric heat of deformation and surface frictional heat flux at material in contact with pin according to Eq. ( 8).In practice, if localized melting happens, the heat generation stops and the extra heat will be absorbed by the bulk material and the temperature decreases below the solidus temperature again.This would result in unwanted weld defects.
As reported in [18,49], there is an empirical softening regime at which the dynamic viscosity becomes almost zero when the temperature reaches near the solidus/ cut-off temperature.We consider a linear decrease of flow stress to zero as temperature increases under different strain rates in the neighboring region of 50 o C below the solidus temperature of aluminum 6061 (582 o C) as shown in Fig. 3.

MESH SENSITIVITY ANALYSIS AND SOLVER SELECTION
The optimum mesh size in the plate shown in elements were used throughout the analysis and as shown in Table 1 and Fig 4, the error (calculated minimum velocity) remains almost constant and near zero when a 'fine' mesh is used.Hence, in order to have the minimum error and the lowest computational cost at the same time, we chose a mesh size between 'fine' and 'finer' in all subsequent cases.Also, in order to have more mesh size control around the pin we divided the part domain into two subdomains as shown in Fig 5 .In case studies 1 to 3 (details to follow in the next section) where we have a pin radius of 5 mm, the total number of meshes were 1758 with 116 boundary elements; and in case 4 where the pin radius is 2.22 mm we used the same number of boundary elements but smaller mesh size around the pin subdomain with a total number of 3132 elements.The chosen solver in COMSOL was UMFPACK with an error tolerance of 10 - 6 which yielded the best results.UMFPACK is a set of routines for solving nonsymmetric sparse linear systems, Ax=b, using the Unsymmetric Multi Frontal method [51].

CASE STUDIES
In order to test the predictability of the developed model against temperature, strain, and strain rate fields via other published reports as well as our temperature measurements during FSW of aluminum 6061 alloy, four different cases with different process conditions were considered as follows.

CASE 1: LONG ET AL [52]: VALIDATION OF THE TEMPERATURE FIELD
In this case, we try to reproduce the results of Long et al research [52] with the same process parameters and material properties as summarized in Table 2.In their 2D steady-state Eulerian model, heat generation was due to viscous dissipation in the material flow (fluid).Viscous dissipation and hence temperature increase was limited by a modification to the viscosity function.The viscosity is forced to decrease exponentially by three orders of magnitude within 50 o C above the solidus temperature of AA6061 (Ts~582 o C), i.e., simulating a steep reduction in the viscosity in the transitional regime from the solid to the semisolid state.This, in turn, prevents the model from attaining temperatures greater than the melting temperature [53].A rigid-viscoplastic material behavior was used according to [54].
When this case was remodeled with our integrated multiphysics model, as explained in the previous section, a linear decrease of flow stress in Eq. ( 10) to zero is considered as the temperature decreases by 50 o C below the solidus temperature, which is slightly different from what was considered in the original model.

CASE 2: EXPERIMENTAL TEMPERATURE MEASUREMENTS: VALIDATION OF THE TEMPERATURE FIELD
In this case, we performed temperature measurements on two points on the back surface of the plates during FSW of aluminum 6061 plates with 6.5 mm thickness, as shown in Fig. 6.
The process parameters and material properties used in the model are summarized in Table 3.

CASE 3: XU ET AL [44]: VALIDATION OF THE STRAIN FIELD
In this case, we remodeled the work of Xu et al [44].In their 2D and 3D steady-state arbitrary Lagrangian-Eulerian models, temperature values measured from an actual FSW process were used as input for the solid mechanics module.Temperature values in the test were measured at mid thickness of the plate with thermocouples inserted in small holes in the plates.At the tool-plate interface, a frictional contact model was proposed to simulate the mechanical interaction between the tool and the plate.Namely, the contact behavior at the tool-plate interface was taken to be governed by a modified Coulomb friction law.Temperature dependent elastic-plastic material properties were also used.In order to capture the detailed distribution of the equivalent plastic strain, a 2D plane-strain model of the horizontal cross-section had also been developed for the same in-plane model geometry as the 3D model.They concluded that, in general, a 2D plane strain model can provide a good approximation of the mechanical fields on the horizontal cross-section at half plate thickness.The same observation was noted when measuring a low strain in thickness direction compared to a high strain in streamlines direction via a combined visioplasticity-CAD approach [43].Here we use the same process parameters (Table 4) as in [44] and the rigidviscoplastic material constants are employed according to Tello et al [50].This material behavior would be more accurate compared to the one used by Sheppard et al [54] and would be closer to the actual hot deformation flow curves as described in [50].

CASE 4: WANG ET AL [55]: VALIDATION OF THE SHEAR STAIN RATE VS. TEMPERATURE
In this case we remodel the temperature and shear strain rate distribution of an Eulerian 3D model at mid thickness of the plate at a central streamline around the pin.In the original work of Wang et al [55], a stick boundary condition was assumed between the tool and the matrix (i.e., __=1).They used the rigid-viscoplastic material model according to Eq.( 10) with the same empirical softening regime explained in section 2.4, in order to calculate the dynamic viscosity of the non-Newtonian hot deformed material.Once the flow rate has been calculated, the heat flux across the pin and shoulder surfaces is found according to Eq (14).
Total heat flux [55]: (14) which is a product of the shear stress at the surface and the wall velocity.The Eq. ( 14) is an empirical relation which was used to calculate the total heat flux on the tool.
Obtaining an accurate estimate of the temperature-dependent thermal conductivity coefficient is difficult as it can be sensitive to the state of the material, i.e. the amount of solute.Handbook values often give the thermal conductivity at equilibrium and therefore are not representative of the thermal conductivity during welding where the thermal cycle is very short.However, when a small volume of material is affected by the FSW process, as indicated by the heat affected zone measurements, the use of room temperature conductivity value may be reasonable [49].The values of the process parameters and the material properties used for this case are given in Table 5. 3.55 [54] 2.41e8 [54] 145000 [54] 789.9+0.4959T[52] 115.23+0.1594T[52] 2700 [52] 582 [52] Table 3: Process parameters and material properties used in Case 2 3.55 [54] 2.41e8 [54] 145000 [54] 900 [3] 160 [3] 2700 [3] 582 [3] Table 4: Process parameters and material properties used in Case 3 5.33 [50] 1.63e13 [50] 191000 [50] 900 [3] 160 [3] 2700 [3] 582 [3] Table 5: Process parameters and material properties used in Case 4

Rotational
Weld speed- 4.70929 [55] 3.0197e11 [55] 168000 [55] 900 [3] 160 [3] 2700 [3] 570 [55] Figure 3: Constitutive behavior of aluminum 6061 based on Eq. ( 10) with constants taken from [50] Figure 4: Effect of mesh size on minimum velocity value (m/s) which is defined as error during mesh sensitivity analysis The result of temperature distribution around the pin in Case 1 is shown in Fig. 7 and compared with the original data from [52].The two temperature fields show a similar distribution in Fig. 7, while the minor difference in maximum temperature is due to different cut-off temperature values assumed in the two models.It should also be noted that in the present model ( As seen in Fig. 8, the maximum temperature is located at the backside of the advancing side of the pin after the material passes the pin and moves along the retreating side.As shown in streamlines of Fig. 8, a shear layer exists around the pin which was proposed by Schneider et al [56].

Rotational
We also examined the effect of using constant (temperature independent) thermal properties in the model and we noted that it generates a higher maximum temperature in the model.Also using the values of the thermo-mechanical model constants , , and n (Eq (10)) from the Sheppard et al study [54] generated a lower maximum temperature when compared to that of Tello et al [50].

CASE 2: EXPERIMENTAL TEMPERATURE MEASUREMENTS
The results of comparison of temperature distributions at points 1 and 2 (Fig. 6) are shown in Fig. 9 via the steady state model.They are plotted using the temperature changes with distance from the tool center, by considering the weld speed and converting temperature-time measurements to temperature-distance values.The model predicts a higher value of the peak temperature (about 30 o C difference compared to the measurement) and also a lower value in areas near the peak (~ 50 o C difference) which could be resulted due to the 2D modeling which does not consider heat generation from the shoulder.

CASE 3: XU ET AL [44]
The result of strain distribution in Case 3 is shown in Fig. 10-a and compared with the original results of Xu et al [44] (Fig 10-b) and a metalworking model from Arbegast [57] (Fig 10-c).The strain distribution at mid thickness of the plate in the present model (Fig. 10a) has almost the same strain distribution as in Fig 10-b.The maximum strain occurs in front of the leading edge of the pin which is estimated at 31.98 by our model.Near the pin there is almost the same/uniform high strain distribution.The maximum strain at mid thickness of  [57].Our model can also predict the formation of onion rings with a space equal to advance per rotation as shown in Fig 10-a.There is a negative strain distribution around the pin in our model which complies with the forging zone in Fig. 10-c.It also indicates the ring vortex generation around the pin as proposed by Schneider et al [56].If the material is trapped in the shear layer, its deposit is delayed and may rotate with the pin more than one cycle.This can explain some of material flow behavior during FSW which has been reported by Colligan [58].4.4 CASE 4: WANG ET AL [55] In case 4, the values of temperature (K) and shear strain rate (s -1 ) are predicted in different regions of a central streamline around the pin (as marked in Fig. 11).The original results from the research of Wang et al [55] with different weld speeds and a constant rotational speed are shown in Fig 12 .The results of our model with the specified process parameters (weld speed of 13.33 mm/sec and the rotational speed of 41.66 rev/sec) are shown in Fig. 13 as well as the results from the research of Wang et al using similar process parameters.Comparing the results, similar patterns and shear strain rate values are seen, while a temperature difference of about 20 to 40 o C in different points exists, which should be the consequence of using different heat models in the two models.Our model predicts slightly less heat generation compared to the work of Wang et al [55] or maybe a slightly different thermal conductivity value is used for the plate and tool which in turn yields slightly lower temperature distribution at the leading edge of the tool (i.e., region 1 in Fig. 11).Figure 14 shows the shear strain rate distribution around the pin.The location of maximum shear strain rate is near the top of the advancing side of the pin.

CONCLUDING REMARKS
The developed integrated multiphysics model coupled different non-Newtonian incompressible fluid dynamics, heat conduction and convention, and plain stress solid mechanics governing relations to model FSW of aluminum 6061.It predicted strain at all points of the work piece for the first time, compared to other CFD-based models which predicted strain only on limited streamlines.The model was successfully applied to predict the distribution of microstructure and residual stress around the pin [40,41], as it predicts Figure 13: Values of temperature (K) and shear strain rate (s-1) in different regions using the present model with the weld speed of 13.33 mm/sec and the rotational speed of 41.66 rev/sec, along with corresponding results from [55] (Fig. 12) temperature distribution at mid thickness with an acceptable tolerance.It can also predict the distribution of temperature, strain, and strain rate around the pin during welding as it was shown through four case studies as well in-house experiments on aluminum 6061.Combining the 'integrated' and "multiphyscis' features under the same model is believed to be the most effective coupled modeling approach for FSW which can be used in future studies to predict interactions of main and post-process parameters for different aluminum alloys.
Based on specific results in each presented case study one may conclude that: 1.
The maximum temperature during FSW is located at the backside of the advancing side of the pin after the material pass as the pin and moves along the retreating side.

2.
The maximum strain during FSW in mid thickness of the plate occurs in front of the leading edge of the pin.

3.
The generation of shear zone, forging zone and ring vortex around the pin can be well explained based on the numerical model results.

4.
The pattern of temperature vs. shear strain rate is predictable in different regions of a central streamline around the pin.The maximum shear strain rate location is near the top of the advancing side of the pin.

5.
If heat transfer in the tool is not considered in the numerical model, a lower value of temperature in the leading edge of the tool is resulted.6.
If constant thermal properties are used in the model, they generate higher maximum temperature.Also using the values of a thermo-mechanical model constants (, , and n) from the study of Sheppard et al [54] generated a lower maximum temperature compared to using those from the work of Tello et al [50].This indicates a necessity for further experimental-numerical studies in the field to arrive at a set of optimum model parameters that would equally perform under different processing conditions.

Figure 1 :
Figure 1: Schematic view of the model dimensions and boundary conditions

Figure 2 :
Figure 2: Schematic view of the weld speed and rotational velocity at the pin surface

= 3 Figure 2 :
Figure 2: Schematic view of the weld speed and rotational velocity at the pin surface

Table 2 :
Process parameters and material properties used in Case 1

Figure 5 :
Figure 5: Mesh size distribution and subdomains Fig 7-a) the tool moves from right to left and rotated counterclockwise but in the work of Long et al (Fig 7-b) [52] the tool moved from left to right and rotated clockwise, though this should not have any effect in the results.Using our multiphysics model, the velocity field streamlines along with the temperature distribution are also shown in Fig. 8.

Figure 6 :Figure 7 :
Figure 6: Points where temperature was measured in the back of the aluminum plates during FSW

Figure 8 :Figure 9 :
Figure 8: Temperature distribution (K) along with the velocity field streamlines in Case 1 (using the present model)

Figure 11 :
Figure 11: Schematic a central streamline with four different regions around the pin in Case 4

Figure 12 :
Figure12: Variation of the temperature and shear strain rate using different welding speeds and a constant rotational speed of 41.66 rev/sec[55]

Figure 13 :
Figure 13: Values of temperature (K) and shear strain rate (s-1) in different regions using the present model with the weld speed of 13.33 mm/sec and the rotational speed of 41.66 rev/sec, along with the corresponding results from [55] (Fig. 12)

Table 1 :
Results of the mesh sensitivity analysis in the model with a pin radius of 5 mm