N. Dialami, M. Chiumenti, M. Cervera and C. Agelet de Saracibar
International Center for Numerical Methods in Engineering (CIMNE)
Universidad Politécnica de Cataluña
c/ Gran Capitán s/n, Módulo C1, Campus Norte UPC, 08034 Barcelona, Spain
This paper deals with the numerical simulation of Friction Stir Welding (FSW) processes. FSW techniques are used in many industrial applications and particularly in the aeronautic and aerospace industries, where the quality of the joining is of essential importance. The analysis is focused either at global level, considering the full component to be jointed, or locally, studying more in detail the heat affected zone (HAZ).
The analysis at global (structural component) level is performed defining the problem in the Lagrangian setting while, at local level, an apropos kinematic framework which makes use of an efficient combination of Lagrangian (pin), Eulerian (metal sheet) and ALE (stirring zone) descriptions for the different computational subdomains is introduced for the numerical modeling. As a result, the analysis can deal with complex (noncylindrical) pinshapes and the extremely large deformation of the material at the HAZ without requiring any remeshing or remapping tools.
A fully coupled thermomechanical framework is proposed for the computational modeling of the FSW processes proposed both at local and global level. A staggered algorithm based on an isothermal fractional step method is introduced.
To account for the isochoric behavior of the material when the temperature range is close to the melting point or due to the predominant deviatoric deformations induced by the viscoplastic response, a mixed finite element technology is introduced. The Variational Multi Scale (VMS) method is used to circumvent the LBB stability condition allowing the use of linear/linear P1/P1 interpolations for displacement (or velocity, ALE/Eulerian formulation) and pressure fields, respectively. The same stabilization strategy is adopted to tackle the instabilities of the temperature field, inherent characteristic of convective dominated problems (thermal analysis in ALE/Eulerian kinematic framework).
At global level, the material behavior is characterized by a thermoelasto viscoplastic constitutive model. The analysis at local level is characterized by a rigid thermoviscoplastic constitutive model. Different thermally coupled (nonNewtonian) fluidlike models as NortonHoff, Carreau or SheppardWright, among others are tested.
To better understand the material flow pattern in the stirring zone, a (Lagrangian based) particle tracing is carried out while postprocessing FSW results.
A coupling strategy between the analysis of the process zone nearby the pintool (local level analysis) and the simulation carried out for the entire structure to be welded (global level analysis) is implemented to accurately predict the temperature histories and, thereby, the residual stresses in FSW.
FSW process. 
Figure 1: FSW process. 
During FSW, the workpiece is placed on a backup plate and it is clamped rigidly to eliminate any degrees of freedom (Figure 1). A nonconsumable tool, rotating at a constant speed, is inserted into the welding line between two pieces of sheet or plate material (which are butted together) and generates heat. This heat is produced, on one hand, by the friction between the tool shoulder and the workpieces, and, on the other hand, by the mechanical mixing (stirring) process in the solid state. This results in the plastification of the material close to the tool at very high strain rates and leads to the formation of the joint. In detail, the plasticized material is stirred by the tool and the heated material is forced to flow around the pin tool to its backside thus filling the hole in the tool wake as the tool moves forward. As the material cools down, a solid continuous joint between the two plates emerges.
Usually the tool is tilted at an angle of away from the direction of travel, although some tool designs allow it to be positioned orthogonally to the workpiece (Figure 2). The welding tool consists of a shoulder and a pin. The length of the pin tool is slightly less than the depth of weld and the tool shoulder is kept in close contact with the workpiece surface (see Figure 3). The tool serves three primary functions, that is, heating of the workpiece, movement of material to produce the joint (stirring), and containment of the hot metal beneath the tool shoulder. The function of the pin tool is to heat up the weld metal by means of friction and plastic dissipation, and, through its shape and rotation, force the metal to move around its form and create a weld. The function of the shoulder is to heat up the metal through friction and to prevent it from being forced out of the weld. The tool shoulder restricts metal flow to a level equivalent to the shoulder position, that is, approximately to the initial workpiece top surface.Tilt angle in FSW 
Figure 2: Tilt angle in FSW 
Error creating thumbnail: File missing

Figure 3: Schematic representation of the friction stir welding process 
Definition of Friction Stir Welding zones Dong00 
Figure 4: Definition of Friction Stir Welding zones Dong00 
Since FSW is not a symmetric process, two sides of the tool are differentiated. One can see in Figure 4 that the workpieces being joined by the weld are either on the retreating or advancing side of the rotating tool. The retreating side is the one where the tool rotating direction is opposite to the tool moving direction and parallel to the metal flow direction. In contrast, the advancing side is the one where the tool rotation direction is the same as the tool moving direction and opposite to the metal flow direction. This unsymmetric nature results in a different material flow on the different sides of the tool and has a large effect on many applications, especially lap joints [27]. The periodic "onion flow" pattern that is left behind as the tool advances is schematically illustrated in Figure 4.
During the early development of FSW, the process appeared simple, compared to many conventional welding practices. However, as development continued, the complexity of FSW was realized. It is now known that properties following FSW are a function of both controlled and uncontrolled variables (response variables) as well as external boundary conditions. For example, investigators have now illustrated that postweld properties depend on:
Different pin shapes. 
Figure 5: Different pin shapes. 
These parameters must be carefully calibrated according to the welding process and the selected material, respectively. The strong coupling between the temperature field and the mechanical behavior is the keypoint in FSW and its highly nonlinear relationship makes the process setup complex. The operative range for most of the welding process parameters is rather narrow requiring a tedious characterization and sensitivity analysis. This is why, despite the apparent simplicity of this novel welding procedure, computational modeling is considered a very helpful tool to understand the leading mechanisms that govern the material behavior, attracting more and more the research interest.
Finite element modeling is an option which can help to determine process parameters that would otherwise require further experimental testing for validation and analysis. The weld quality depends largely on how the material is heated, cooled and deformed. Hence a prior knowledge of the temperature evolution within the workpiece would help in designing the process parameters for a welding application. Research in the field of FSW joints has been limited possibly due to proprietary publishing restriction within industry. For this reason, Finite Element Analysis (FEA) could be also very beneficial. Two process parameters of interest for FSW welds are tool travel rates and rotational tool velocities. With respect to this, a lot of emphasis has been laid on FEA analysis, as it may broaden the scope of application of FSW. Another important process parameter in FSW is the heat flux. This can be also easily included in the FEA. The heat flux should be high enough to keep the maximum temperature in the workpiece around 80% to 90% of the melting temperature of the workpiece material, so that welding defects are avoided [33].
Moreover, analytical and numerical methods have a role to play, although numerical methods dominate due to the accuracy and easetouse of modern workstations and software. Numerical modeling is based on discretized representations of specific welds, using finite element, finite difference, or finite volume techniques. These methods can capture much of the complexity in material constitutive behavior, boundary conditions, and geometry, but in practice, a limited range of conditions tends to be studied in depth. Therefore, it is good modeling practice to explore simplifications to the problem that give useful insight across a wider domain, for example, making valid twodimensional (2D) approximations to inherently three dimensional (3D) behavior. It is also essential to deliver a model that is properly validated and whose sensitivity is known to uncertainty in the input material and process data–ideals that are rarely carried through in practice.
(^{1}) However, it must be noted that, as in the traditional fusion welds, there also exist a softened heat affected zone (HAZ) and a tensile residual stress parallel to the weld.
Information about the shape, dimensions and residual stresses in a component after welding and mechanical properties of the welded joint are of great interest in order to improve the quality and to prevent failures during manufacturing or in service. The FSW process can be analyzed either experimentally or numerically.
FSW is difficult to analyze experimentally; however, process parameters and different fixture setups can be evaluated without doing a large number of experiments. An experiment can be designed to answer one or more carefully formulated questions. The goals must be clarified perfectly to choose the appropriate parameters and factors. Otherwise the goal is not achieved and the experiment must be repeated. Different welded specimens are produced by varying the process parameter. The properties and microstructure changes in weld are investigated. For instance, the tensile strength of the produced joint is tested at room temperature. Microstructure of the weld is analyzed by means of optical microscopy or microhardness measurements.
The alternative to the experimental FSW analysis is numerical modeling and simulation. Computerbased models provide the opportunity to improve theories of design and increase their acceptance. Simulations are useful in designing the manufacturing process as well as the manufactured component itself. To do an appropriate modeling, the physics of the problem must be wellunderstood.
FSW simulation is a problem of complex nature; the process is highly nonlinear and coupled. Different physical phenomena occur during the welding process, involving the thermal and mechanical interactions. The temperature field is a function of many welding parameters such as welding speed, welding sequence and environmental conditions. Formation of distortions and residual stresses in workpieces depend on many interrelated factors such as thermal field, material properties, structural boundary conditions and welding conditions.
The challenging issues in physical modeling of the FSW process are divided into three parts:
Heat transfer mechanisms including convection, radiation and conduction have a significant role on the process behavior. Convection and radiation fluxes dissipate heat significantly through the workpieces to the surrounding environment while conduction heat flux occurs between the workpieces and the support.
In FSW, the mechanical behavior is nonlinear due to the high strain rates and viscoplastic material. The strong nonlinear region is limited to a small area and the remaining part of the model is mostly linear. However, the exact boundaries of the nonlinear zone are not known apriori. It is generally believed that strain rate during the welding is high. Knowledge of strain rate is important for understanding the subsequent evolution of grain structure, and it serves as a basis for verification of various models as well.
The thermal and mechanical problems are strongly coupled (the thermal loads generate changes in the mechanical field). The mechanical effects coupled to the thermal ones include internal heat generation due to plastic deformations or viscous effects, heat transfer between contacting bodies, heat generation due to friction, etc. The thermal effects are also coupled to the mechanical ones; for instance, thermal expansion, temperaturedependent mechanical properties, temperature gradients in workpieces, etc.
An adequate physical model of the welding process must account for all these phenomena including thermal, mechanical and coupling aspects.
Among several numerical modeling techniques, the Finite Element (FE) framework is found to be suitable for the simulation of welding and proven to be a versatile tool for predicting a component's response to the various thermal and mechanical loads. The FE method also offers the possibility to examine different aspects of the manufacturing process without having a physical prototype of the product. To this end, a specialized thermomechanical coupled model needs to be implemented in a finite element program, and the predictive capabilities of the theory and its numerical implementation must be validated.
The numerical simulation of the FSW process has many complex and challenging aspects that are difficult to deal with: the welding process is described by the equilibrium and energy equations governing the mechanical and thermal problems and they are coupled. Additionally, both of them are nonlinear. This has important implications upon the complexity of the numerical model. Consequently, a robust and efficient numerical strategy is crucial for solving such highly nonlinear coupled finite element equations. In such process, several assumptions are commonly assumed.
It is important to distinguish between two different kinds of welding analyses carried out at local or global level, respectively.
In local level analysis, the focus of the simulation is the heat affected zone. The simulation is intended to compute the heat power generated either by viscoplastic dissipation or by friction at the contact interface. At this level, the process phenomena that can be studied are the relationship between welding parameters, the contact mechanisms in terms of applied normal pressure and friction coefficient, the setting geometry, the material flow within the heat affected zone, its size and the corresponding consequences on the microstructure evolution, etc.
A simulation carried out at global level studies the entire component to be welded. In this case, a moving heat power source is applied to a control volume representing the actual heat affected zone at each timestep of the analysis. The effects induced by the welding process on the structural behavior are the target of this kind of study. These effects are distortions, residual stresses or weaknesses along the welding line, among others.
The aim of this work is to develop a robust numerical tool able to simulate the welding process considering its complex features at local level as well as global level.
A quasistatic mechanical analysis can be assumed as the inertia effects in welding processes are negligible due to the high viscosity characterization. At local level, the volumetric changes are found to be negligible, and incompressibility can be assumed. To deal with the incompressible behavior, a very convenient and common choice is to describe the formulation splitting the stress tensor into its deviatoric and volumetric parts. Dealing with the incompressible limit requires the use of mixed velocitypressure interpolations. The problem suffers from instability if the standard Galerkin FE formulation is used, unless compatible spaces for the pressure and the velocity field are selected (LadyzhenskayaBabuska–Brezzi (LBB) stability condition). Due to this, pressure instabilities appear if equal velocitypressure interpolations are used. Thus, the challenging issue of pressure stabilization rises up.
The welding process is characterized by very high strain rates as well as a wide temperature range going from the environmental temperature to the melting point. Hence, the constitutive laws adopted should depend on both variables. The constitutive theory applied must be specialized to capture the features of the thermomechanically coupled strain rate and temperature dependent large deformation. According to the split of the stress tensor, different ratedependent constitutive models can be used for modeling of the welding process. At typical welding temperatures, the large strain deformation is mainly viscoplastic. Depending on the scope of the analysis, rigidviscoplastic or elastoviscoplastic constitutive models can be used. Not only the prediction of the temperature evolution, but the accurate residual stress evaluation field generated during the process is the objective of the FSW simulation. The selected constitutive model must appropriately define the material behavior and has to be calibrated by the temperature evolution. The challenge arises from the extremely nonlinear behavior of these constitutive models and, therefore, from the numerical point of view, a special treatment is obligatory. Moreover, the localized large strain rates usually involved in FSW processes make the problem even more complex.
The thermal problem is defined by the balance of energy equation. In FSW simulation, the plastic dissipation term appearing in the energy equation has a critical role on the process behavior and it is the main source of heat generation.
The definition of the heat source is one of the key points when studying the welding process. In global level simulations, the mesh density used to discretize the geometry is not usually fine enough to define the welding pool shape or a nonuniform heat source. This is only done if the simulation of the welding pool is the objective itself (local level analysis). If the global structure is considered, the size of the heat source is of the same dimension than the element size generally used for a thermomechanical analysis. Therefore, when the global model is taken into account, the resulting mesh density is usually too coarse to represent the actual shape of the heat source.
Depending on the framework used to describe the formulation of the coupled thermomechanical problem, a convective term might appear in the thermal governing equations. Therefore convection instabilities of the temperature appear for convection dominated problems. It is well known that in diffusion dominated problems, the solution is stable. However, in convection dominated problems, the stabilizing effect of the diffusion term becomes insufficient and oscillations appear in the temperature field. The threshold between stable and unstable solutions is usually expressed in terms of the Peclet number.
Establishing an appropriate kinematic framework for the simulation of welding is one of the main objectives of this paper.
If the welding process is studied at global level, a Lagrangian framework is an appropriate choice for the description of the problem. Lagrangian settings, in which each individual node of the computational mesh represents an associated material particle during motion, are mainly used in structural mechanics. Classical applications of the Lagrangian description in large deformation problems are the simulation of vehicle crash tests and the modeling of metal forming operations. In these processes, the Lagrangian description is used in combination with both solid and structural (beam, plate, shell) elements. Numerical solutions are often characterized by large displacements and deformations and historydependent constitutive relations are employed to describe elastoplastic and viscoplastic material behavior. The Lagrangian reference frame allows easy tracking of free surfaces and interfaces between different materials.
In a local simulation, the main focus of the simulation is the Heat Affected Zone (HAZ) where the use of a Lagrangian framework is not always advantageous. In the HAZ, the large distortions would require continuous remeshing. The alternative is to use Eulerian or Arbitrary Lagrangian Eulerian (ALE) methods. Eulerian settings are widely used in fluid mechanics. Here, the computational domain and reference mesh are fixed and the fluid moves with respect to the grid. The Eulerian formulation facilitates the treatment of large distortions in the fluid motion. Its handicap is the difficulty to follow free surfaces and interfaces between different materials or different media (e.g., fluidfluid and fluidsolid interfaces).
An Arbitrary Lagrangian Eulerian (ALE) formulation which generalizes the classical Lagrangian and Eulerian descriptions is particularly useful in flow problems involving large distortions in the presence of mobile and deforming boundaries. Typical examples are problems describing the interaction between a fluid and a flexible structure and the simulation of metal forming processes. The key idea in the ALE formulation is the introduction of a computational mesh which can move with a velocity different from (but related to) the velocity of the material particles. With this additional freedom with respect to the Eulerian and Lagrangian descriptions, the ALE method succeeds to a certain extend in minimizing the problems encountered in the classical kinematical descriptions, while combining their respective advantages at best.
In the simulation of FSW, it is adroit to introduce an apropos kinematic framework for the description of different parts of the computational domain. Despite the efficiency of the idea, the mesh moving strategy and the treatment of the domains interaction are challenging.
The numerical solution of the coupled thermomechanical problem involves the transformation of an infinite dimensional transient system into a sequence of discrete nonlinear algebraic problems. This is achieved by means of the FE spatial descritization procedure, a timemarching scheme for the advancement of the primary nodal variables, and with a time iteration algorithm to update the internal variables of the constitutive equations.
Regarding the timestepping schemes, two types of strategies can be applied to the solution of the coupled thermomechanical problems:
The first possibility is a monolithic (simultaneous) timestepping algorithm which solves both the mechanical and the thermal equilibrium equations together. It advances all the primary nodal variables of the problem simultaneously. The main advantage of this method is that it enables stability and convergence of the whole coupled problem. However, in simultaneous solution procedures, the timestep as well as timestepping algorithm has to be equal for all subproblems, which may be inefficient if different time scales are involved in the thermal and the mechanical problem. Another important disadvantage is the considerably high computational effort required to solve the monolithic algebraic system and the necessity to develop software and solution methods specifically for each coupled problem.
A second possibility is a staggered algorithm (blockiterative or fractionalstep), where the two subproblems are solved sequentially. Usually, a staggered solution (arising from an operator split and a product formula algorithm (PFA)) yields superior computational efficiency.
Staggered solutions are based on an operator split, applied to the coupled system of nonlinear ordinary differential equations, and a product formula algorithm, which leads to splitting of the original monolithic problem into two smaller and better conditioned subproblems (within the framework of classical fractional step methods). This leads to the partition of the original problem into smaller and typically symmetric (physical) subproblems. After this, the use of different standard timestepping algorithms developed for the uncoupled subproblems is straightforward, and it is possible to take advantage of the different time scales involved. The major drawback of these methods is the possible loss of accuracy and stability. However, it is possible to obtain unconditionally stable schemes using this approach, providing that the operator split preserves the underlying dissipative structure of the original problem.
One of the main issues in the study of FSW at local level, is the heat generation. The generated heat must be enough to allow for the material to flow and to obtain a deep heat affected zone. Insufficient heat forms voids as the material is not softened enough to flow properly. The visualization of the material flow is very useful to understand its behavior during the weld. A method approving the quality of the created weld by visualization of the joint pattern is advantageous. It can be used to investigate the appropriate process parameters to create a qualified joint. However, following the position of the material during the welding process is not an easy task, neither experimentally or numerically.
The experimental material visualization is difficult and needs metallographic tools. This is why establishing a numerical method for the visualization of the material trajectory in order to gain insight to the heat affected zone and the material penetration within the workpiece thickness is one of the main objective of the work. Particle tracing is a method used to simulate the motion of material points, following their positions at each timestep of the analysis. This method can be naturally applied to the study of the material flow in the welding process. In the Lagrangian framework, as the mesh nodes represent the material points, the trajectories are the solution of the governing system of equations. When using Eulerian and ALE framework the solution does not give directly information about the material position. However, the velocity field obtained can be integrated to get an insight of the extent of material mixing during the weld.
Integration of the velocity field is proposed at postprocess level to follow the material motion (displacement field). This obliges the modeler to use an appropriate time integration method for the solution of the ODE in order to track the particles. Moreover, a search algorithm must be executed to find the position of the material points in the Eulerian or ALE meshes.
Since FSW occurs by the deformation of material at temperatures below the melting point, many problems commonly associated with fusion welding technologies can be avoided and highquality welds are produced.
Generally, FSW yields fine microstructures, absence of cracking, low residual distortion and no loss of alloying elements. Nevertheless, as in the traditional fusion welds, a softened heat affected zone and a tensile residual stress field appear.
Although the residual stresses and distortion are smaller in comparison with those of traditional fusion welding, they cannot be ignored, specially when welding thin plates of large size.
In the local level analysis, the focus of the study is the HAZ and a viscoplastic model is used to chareacterize the material behavior. Elastic stresses are neglected and the calculation of residual stresses is not possible. However, at global level, the residual stresses are one of the main outcomes of the process simulation using an elastoviscoplastic constitutive model. Therefore, in this work, a localglobal coupling strategy is proposed in order to obtain the residual stress field.
In FSW the tool generates heat via a combination of friction at the tool surface and viscous dissipation within the deformed material. This heat is conducted into the tool and welded material, and is then convected from the top surface or conducted into the backing plate. The majority of both analytical and numerical models are used to describe this heat flow and to see its effect on the entire piece during and after FSW. The following section lists some works performed to model FSW. They are divided in: analytical and numerical modeling of FSW, thermal and thermomehanical modeling, global analysis and local modeling of the heat input and use of different kinematic frameworks.
Thermal models are used to describe the heat flow in the FSW process. Midling and Grong [83] presented a transient analytical model that used a modified version of Rykalin's infinite rod solution and predicted the weld temperature and microstructure. The microstructural evolution was predicted for 6082T6 and AISiC, and indicated that a narrow HAZ was obtained when a high specific power was used in conjunction with a short duration heating cycle.
Grong [66] described how Rosenthal's [98] analytical solution to the heat equation could be used to predict the thermal profile. Three different analytical solutions existed for thin, medium and thick plate. The appropriate solution was determined by the power and speed of the heat source, and the physical dimensions and thermal properties of the material being welded using a dimensionless heat flow map. When using the medium plate solution, it was necessary to superimpose virtual heat sources to include the effect of plate thickness and width. Finally, analytical solutions assumed constant thermal properties.
Early work by Russell and Shercliff [99], [Russell Sheercliff 1999] used a point heat source which was considered adequate because the temperature in the HAZ was of interest. The analytical solution was extended by McClure et al. [81], and Gould and Feng [65] by distributing the heat source over the shoulder. The distributed heat source was integrated to find the temperature at a particular position in the material. The heat source intensity increased with the distance from the centre of the tool in McClure et al. [81] and was imposed on a ring around the shoulder in Gould and Feng [65]. Imposing the heat source as a ring resulted in a large temperature peak around the shoulder. When the heat input was correctly adjusted, reasonable correlation was obtained with experimental results.
Finally, an inverse problem approach was used by Fonda and Lambrakos [58]. Rather than using thermocouple measurements, they examined the weld microstructure and assumed that the edge of the heat affected zone corresponded to a peak weld temperature of 250 C. Using this knowledge, they worked backwards to find the heat input and weld thermal profile at all points in the weld. Stewart et al. [112] also used an analytic thermal model, however the main emphasis of this work was material flow.
Song and Kovacevic [111] developed a mathematical model to describe the detailed threedimensional transient heat transfer process. Their work was both theoretical and experimental. An explicit central differential scheme was used in solving the control equations. The heat input from the tool shoulder was modeled as frictional heat and the heat from the tool pin as uniform volumetric heat generated by the plastic deformation near the pin. Moving coordinates and a nonuniform grid mesh were introduced to reduce the difficulty of modeling the heat generation due to the movement of the tool pin.
In the recent work of Ferrer et. al [57], a simple analytical solution using series expansions for Cauchy momentum and energy equationset is obtained. The Power Law model takes into account the shear thinning fluid and the Arrheniustype relationship the temperature dependent viscosity. Friction dissipation as an external heat source is considered.
As previously discussed, numerical modeling and simulation is an important tool for understanding the mechanisms of the FSW process. It enables obtaining both the qualitative and quantitative insight of the welding characteristics without performing costly experiments. Flexibility of the numerical methods (in particular FEM) in treating complex geometries and boundary conditions defines an important advantage of these methods against their analytical counterparts.
The FSW simulation typically involves studies of the transient temperature and its dependence on the rotation and advancing speed, residual stresses in the workpiece, etc. This simulation is not an easy task since it involves the interaction of thermal, mechanical and metallurgical phenomena. Up to now several researchers have carried out computational modeling of FSW.
To date, most of the research interest devoted to the topic was focused on the heat transfer and thermal analysis in FSW while the mechanical aspects were neglected. Among others, Gould and Feng [65] proposed a simple heat transfer model to predict the temperature distribution in the workpiece. Chao and Qi [31], [32] developed a moving heat source model in a finite element analysis and simulated the transient evolution of the temperature field, residual stresses and residual distortions induced by the FSW process. Their model was based on the assumption that the heat generation came from sliding friction between the tool and material. This was done by using Coulomb's law to estimate the friction force. Moreover, the pressure at the tool surface was set constant and thereby enabled a radially dependent surface heat flux distribution generated by the tool shoulder. In this model the heat generated by the pin was neglected.
Nguyen and Weckman [88] demonstrated a transient thermal FEM model for friction welding which was used to predict the microstructure of 1045 steel. Measured power data was used to calculate the heat input and a constant temperature boundary condition at the welding interface was invoked. In another FEM thermal model by Mitelea and Radu [85] friction welding of dissimilar materials was modeled. The paper compared different heat flux distributions to determine which gave the best agreement with experimental results. To conclude both analytical and numerical techniques were used to describe the heat flow in friction welding, with a more accurate solution being obtained with the latter. Because friction welding is a short duration process, a transient model was more appropriate than one that used a steady state solution.
Colegrove et al. [44], [45] and Frigaard et al. [62] developed 3D heat flow models for the prediction of the temperature field. They used the CFD commercial software FLUENT for a 2D and 3D numerical investigation on the influence of pin geometry, comparing different pin shapes in terms of their influence upon the material flow and welding forces on the basis of both stick and slip conditions at the tool/workpiece interface. It was only the tool pin that was modeled. Several different tool shapes were considered. The modeling result showed that the difference between the result corresponding to slip and stick conditions was small and the pressure and forces were similar. In spite of the good obtained results, the accuracy of the analysis was limited by the assumption of isothermal conditions. Midling [84] and Russell and Sheercliff [100] investigated the effect of tool shoulder of the pin tool on the heat generation during the FSW operation. Generally, those early flow models and others (e.g. Askari et al. [8]) were uncoupled or sequentially coupled to the heat solvers, and limited by the computational power and software capabilities of that time.
More recently, a coupled thermomechanical modeling and simulation of the FSW process can be found in Zhu and Chao [122], Jorge Jr. and Balancín [74] or De Vuyst et al. [49], [48]. Zhu et. al [122] used a 3D nonlinear thermal and thermomechanical numerical model using the finite element analysis code WELDSIM. The objective was to study the variation of transient temperature and residual stress in a friction stir welded plate of 304L stainless steel. Based on the experimental records of transient temperatures, an inverse analysis method for thermal numerical simulation was developed. After the transient temperature field was determined, the residual stresses in the welded plate were calculated using a threedimensional elastoplastic thermomechanical model. In this model the plastic deformation of the material was assumed to follow the Von Mises yield criterion and the associated flow rule.
In a more sophisticated way, De Vuyst et al. [49], [48] used the coupled thermomechanical FE code MORFEO to simulate the flow around tools of simplified geometry. The rotation and advancing speed of the tool were modeled using prescribed velocity fields. An attempt to consider features associated to the geometrical details of the probe and shoulder, which had not been discretized in the FE model in order to avoid very large meshes, was taken into account using additional velocity boundary conditions. In spite of that, the mesh used resulted to be large: a mesh of roughly 250,000 nodes and almost 1.5 million of linear tetrahedral elements was used. A NortonHoff rigidviscoplastic constitutive equation was considered, with averaged values of the consistency and strain rate sensitivity constitutive parameters determined from hot torsion tests performed over a range of temperatures and strain rates. The computed streamlines were compared with the flow visualization experimental results obtained using copper marker material sheets inserted transversally or longitudinally to the weld line. The simulation results correlated well when compared to markers inserted transversely to the welding direction. However, when compared to a marker inserted along the weld center line only qualitative match could be obtained. The correlation could have been improved by modeling the effective weld thickness of the experiment, using a more realistic material model, for instance, by incorporating a yield stress or temperature dependent properties, more exact prescription of the velocity boundary conditions or refining the mesh in specific zones, for instance, under the probe. The authors concluded that it was essential to take into account the effects of the probe thread and shoulder thread in order to get realistic flow fields.
Fourment [59] performed the simulation of transient phases of FSW with FEM using an ALE formulation, in order to take the large deformation of a 3D coupled thermomechanical model into account. This method permitted both transient and steady state analysis. The formulation was developed in [60] and [67] to simulate the different stages of the FSW process. Assidi [10] presented a 3D FSW simulation based on friction models calibration using Eulerian and ALE formulation. An interesting comparison of the heat energy generated by the FSW between numerical methods and experimental data was presented in Dong et al. [55] and Chao et al. [33].
In [33], Chao used a FE formulation to model the heat transfer of the FSW process in two boundary values problems: a steady state problem for the tool and transient one for the workpiece. To validate the result, the temperature evolution was recorded in the tool and in the workpiece. The heat input from the tool shoulder was assumed to be linearly proportional to the distance from the center of the tool due to heat generation by friction. To model the workpiece the code WELDSIM was used. It was a transient, nonlinear, 3D FE code. In this model only half of the workpiece was modeled due to symmetry meaning that the advancing and retreating sides of the weld were not differentiated. The conclusions from this work indicated that 95 % of the heat generated goes into the workpiece and only 5% goes to the tool. It gave a very high heat efficiency estimate.
As the model predictions are not always in agreement with experimental results. In [93], the LevenbergMarquardt (LM) method is used in order to perform a nonlinear estimation of the unknown parameters present in the heat transfer and fluid flow models, by adjusting the temperatures results obtained with the models to temperature experimental measurements. The unknown parameters are: the friction coefficient and the amount of adhesion of material to the surface of the tool, the heat transfer coefficient on the bottom surface and the amount of viscous dissipation converted into heat.
Most of the abovementioned works were performed at global level. These studies typically analyzed the effect of the welding process on the structural behavior in terms of distortions, residual stresses or weakness along the welding line, among others. As the simulations carried out at global level consider the generated heat as input parameter, several techniques were used to determine the heat input to the model. One way was to measure the temperature experimentally, and adjust the heat input of the model till the numerical and experimental temperature profiles match [108], [109], [110], [31], [ChaoQi 1999], [44]. Another technique involved estimating the power input analytically and then artificially limiting the peak temperature [108], [62] or introducing latent heat effects [108], [109] to avoid overpredicting the weld temperature. The most satisfactory approach involved measuring the weld power experimentally and using this as an input to the model [79], [107], [75], [76], [77].
Khandkar et al. [76] used a finite element method based on a 3D thermal model to study the temperature distributions during the FSW process. The moving heat source generated by the rotation and linear traverse of the pintool was correlated to input torque data obtained from experimental investigation of buttwelding. The moving heat source included heat generation due to torques at the interface between the tool shoulder and the workpiece, the horizontal interface between the pin bottom and the workpiece, and the vertical interface between the cylindrical pin surface and the workpiece. Temperaturedependent properties of the weldmaterial were used for the numerical modeling.
In Khandkar et al. [77] and Hamilton et al. [68] a torquebased heat input was used. Various aluminium alloys were included into the model and the maximum welding temperature could be predicted from tool geometry, welding parameters and material parameters. The thermal model involved an energyslip factor which was developed by a relationship between the solidus temperature and the energy per unit length of the weld. In Khandkar et al. latest models, the thermal models were coupled with the mechanical behavior and thereby not only the heat transport was modeled, the residual stress was also an outcome of the model [78].Khandkar et al. [78] used a coupled thermomechanical FE model based on torque input for calculating temperature and residual stresses in aluminium alloys and 304L stainless steel.
Reynolds et al. used two models in [97] to explain the FSW process. The first was a thermal model to simulate temperature profiles in friction stir welds. The total torque at the shoulder was divided into shoulder, pin bottom and vertical pin surface. The required inputs for the model were total input power, tool geometry, thermophysical properties of the material being welded, welding speed and boundary conditions. The output from the model could be used to rationalize observed hardness and microstructure distributions. The second model was a fully coupled, twodimensional fluid dynamics based model that was used to make parametric studies of variations in properties of the material to be welded (mechanical and thermophysical) and variations in welding parameters. This was done by a non slip boundary condition at the tool workpiece interface. The deformation behavior was based on deviatoric flow stress using the ZenerHollomon parameter. Results from this model provided insight regarding the effect of material properties on friction stir weldability and on potential mechanisms of defect formation.
Some other authors presented thermomechanical models for the prediction of the distribution of the residual stresses in the process of friction stir welding. A steadystate simulation of FSW was carried out by Bastier et al. [12]. The simulation included two main steps. The first one uses an Eulerian description of the thermomechanical problem together with a steadystate algorithm detailed in [47], in order to avoid remeshing due to the pin motion. In the second step, a steadystate algorithm based on an elastoviscoplastic constitutive law was used to estimate the residual state induced by the process.
Some other authors used both experimental and numerical methods for computing the residual stresses. McCune et al. [82] studied computationally and experimentally the effect of FSW improvements in terms of panel weight and manufacturing cost on the prediction of residual stress and distortion in order to determine the minimum required modeling fidelity for airframe assembly simulations. They proved the importance of accurately representing the welding forging force and the process speed.
Paulo et al. [92] used a numericalexperimental procedure (contour method) to predict the residual stresses arising from FSW operations on stiffened panels. The contour method allowed for the evaluation of the normal residual stress distribution on a specimen section. The residual stress distribution was evaluated by means of an elastic finite element model of a cut sample, using the measured and digitalized outofplane displacements as input nodal boundary conditions.
Yan et al. [118] adopted a general method with several stiffeners designed on the sheet before welding. Based on the numerical simulation of the process for sheet with stiffeners, the residual distortion of the structure was predicted and the effect of the stiffeners was investigated. They verified first the numerical model experimentally and then applied the verified model on the structure to compute the residual stresses.
Fratini and Pasta [61] used the cutcompliance and the inverse weightfunction methodologies for skin stringer FSW geometries via finite element analysis to measure residual stresses.
Rahmati Darvazi and Iranmanesh [96] presented a thermomechanical model to predict the longitudinal residual stress applying a socalled advancing retreating factor. The uncoupled thermomechanically equations were solved using ABAQUS.
Since the temperature is crucial in the FSW simulation, the heat source needs to be modeled accurately. This consideration obliges the researchers to study the process at the local level where the simulation is concentrated on the stirring zone. For this type of simulation, the heat power is assumed to be generated either by the viscoplastic dissipation or by the friction at the contact interface. At this level the majority of the process phenomena can be analyzed: the relationship between rotation and advancing speed, the contact mechanisms, the effect of pin shape, the material flow within stir zone, the size of the stir zone and the corresponding consequences on the microstructure evolution, etc.
Most models of FSW consider the FEM in a Lagrangian mesh [80], [64] and [63], which is used to study the process globally and to predict the weld temperature and deformation structure. However, the number of simulations using other numerical methods such as finite volume or other kinematic frameworks such as Eulerian or ALE is also considerable.
In the model by Maol and Massoni [80] (FEM with a Lagrangian mesh), the material was assumed to be viscoplastic, with temperaturedependent properties. A frictional interface was used between the two parts and its value was determined experimentally from the pressure and velocity between the two parts. Bendzsak and North [15] used the finite volume method to predict the flow field in the fully plasticized region of a friction weld. Similar and dissimilar welds were analyzed. In the similar welds the viscosity was found from a heuristic relationship, which was independent of temperature. The dissimilar welds used a transient thermal model, a complex viscosity relationship and an EulerianLagrangian mesh.
A thermomechanical model for FSW was proposed by Dong et al. [54]. There, an axissymmetrical FE Lagrangian formulation was used. Ulysse in [113] modeled 3D FSW for aluminium thick plates. Forces acting on the tool were studied for various welding and rotational speeds. The deviatoric stress tensor was used by Ulysse to model the stirwelding process using 3D viscoplastic modeling. Parametric studies were conducted to determine the effect of tool speed on plate temperatures and to validate the model predictions by comparing with available measurements. In addition, forces acting on the tool were computed for various welding and rotational speeds. It was found that pin forces increased with increasing welding speeds, but the opposite effect was observed for increasing rotational speeds.
Askari et al. [9] used the CTH finite volume hydrocode coupled to an advectiondiffusion solver for the energy balance equation. This model predicts important fields like strain, strain rate and temperature distribution. The elastic response was taken into account in this case. The results proved encouraging with respect to gaining an understanding of the material flow around the tool. However, simplified friction conditions were used. They used particle tracking and the mixing fraction to visualize the flow. The mixing fraction determines the ratio of advancing to retreating side material in the welded joint. These techniques enabled very impressive flow visualization, with the large dispersion of material that occurs with an advancing side marker being correctly predicted. The work also predicted that material, particularly that starting on the advancing side, flows more than one revolution around the tool.
Xu et al. [115] used finite element models to describe the material flow around the pin. This was done by using a solid mechanical 2D FE model. It included heat transfer, material flow, and continuum mechanics. The pin was included but not the threads of the pin. Xu and Deng [116], [117] developed a 3D finite element procedure to simulate the FSW process using the commercial FEM code ABAQUS, focusing on the velocity field, the material flow characteristics and the equivalent plastic strain distribution. The authors used an ALE formulation with adaptive meshing and considered large elastoplastic deformations and temperature dependent material properties. However, the authors did not perform a fully coupled thermomechanical simulation, super imposing the temperature map obtained from the experiments as a prescribed temperature field to perform the mechanical analysis. The numerical results were compared to experimental data available, showing a reasonably good correlation between the equivalent plastic strain distributions and the distribution of the microstructure zones in the weld. They examined the velocity gradient around the pin and found that it was higher on the advancing than retreating side and higher in front than behind the pin. Maps of the strain on the transverse cross section were compared against typical weld macrosections.
Seidel and Reynolds [104] also used the CFD commercial software FLUENT to model the 2D steadystate flow around a cylindrical tool. The paper describes the progressive development of a finite volume model in FLUENT that used a viscoplastic material whose viscous properties were based on the SellarsTegart relationship. Even though the model was 2D, heat generation and conduction were included. To avoid overpredicting the weld temperature, the viscosity was reduced by 3 orders of magnitude near the solidus. The model correctly predicted material flow around the retreating side of the tool.
Bendzsak et al. [13], [14] used the Eulerian code Stir3D to model the flow around a FSW tool, including the tool thread and tilt angle in the tool geometry and obtaining complex flow patterns. The temperature effects on the viscosity were neglected. They used the finite volume method to predict the flow field in the fully plasticized region of a friction weld. Similar and dissimilar welds were analyzed. In the similar welds the viscosity was found from a heuristic relationship, which was independent of temperature. The dissimilar welds used a transient thermal model, a complex viscosity relationship and an EulerianLagrangian mesh.
Schmidt and Hattel [103] presented a 3D fully coupled thermomechanical FE model in ABAQUS/Explicit using the ALE formulation and the JohnsonCook material law. The flexibility of the FSW machine was taken into account connecting the rigid tool to a spring. The workpiece was modeled as a cylindrical volume with inlet and outlet boundary conditions. A rigid backplate was used. The contact forces were modeled using a Coulomb friction law, and the surface was allowed to separate. Heat generated by friction and plastic deformation was considered. The simulation modeled the dwell and weld phases of the process. A constant contact conductance was used everywhere under the sheet. They used the generated heat from three different areas, the shoulder, pin and pin tip, and used both sticking and sliding conditions. Despite the wealth of information that this model can provide (e.g. material velocity, plastic strains, and temperatures across the weld), a major shortcoming for it was the long processing time for reaching the steadystate (14 days on a 3 GHz Pentium PC to reach only 10 seconds of model time).
The model developed by Chen and Kovacevic in [34] uses the commercial FEM software ANSYS for a thermomechanically coupled Lagrangian finite element model of aluminium alloy AA6061T6. The welding tool was modeled as a heat source. The model only included the shoulder and so the effect of the pin was ignored. This simple model severely limited the accuracy of the stress and force and the strain rate dependence was not included in the material model. However, the authors were able to investigate the effect of the heat moving source on the workpiece material. Finally the model predicted the welding forces in the x, y and z directions.
Nikiforakis [89] used a finite difference method to model the FSW process. Despite the fact that he was only presenting 2D results, the model proposed had the advantage of minimizing calibration of model parameters, taking into account a maximum of physical effects. A transient and fully coupled thermofluid analysis was performed. The rotation of the tool was handled through the use of the overlapping grid method. A rigidviscoplastic material law was used and sticking contact at the tool workpiece interface was assumed. Hence, heating was due to plastic deformation only.
Heurtier et al. [69] used a 3D semianalytical coupled thermomechanical FE model to simulate FSW processes. The model uses an analytical velocity field and considers heat input from the tool shoulder and plastic strain of the bulk material. Trajectories, temperature, strain, strain rate fields and microhardness in various weld zones are computed and compared to experimental results obtained on a AA 2024T351 alloy FSW joint.
Among the distinctive local level studies, Cho et al. [43] used an Eulerian approach including thermomechanical models without considering the transient temperature in simulation. The strain hardening and texture evolution in the friction stir welds of stainless steel was studied in this paper. A Lagrangian approach with intensive remeshing was employed in [35] while similar approaches were applied in [19] and [20], which are not numerically efficient.
Buffa et al. [19] using the commercial FE software DEFORM3D, proposed a 3D Lagrangian, implicit, coupled thermomechanical numerical model for the simulation of FSW processes, using a rigidviscoplastic material description and a continuum assumption for the weld seam. The proposed model is able to predict the effect of process parameters on process variables, such as the temperature, strain and strain rate fields, as well as material flow and forces. Reasonably good agreement between the numerically predicted results, on forces and temperature distribution, and experimental data was obtained. The authors found that the temperature distribution about the weld line is nearly symmetric because the heat generation during FSW is dominated by rotating speed of the tool, which is much higher than the advancing speed. On the other hand, the material flow in the weld zone is non symmetrically distributed about the weld line because the material flow during FSW is mainly controlled by both advancing and rotating speeds.
Nandan et al. in [86] and [87] employed a control volume approach for discretization of the FSW domain. They investigated viscoplastic flow and heat transfer during friction stir welding in mild steel. The temperature, cooling rates and plastic flows were solved by the equations of conservation of mass, momentum and energy together with the boundary conditions. In this model the nonNewtonian viscosity was determined from the computed values of strain rate, temperature and material properties. Temperatures and total torque was compared with experimental values showing good agreement. The computed temperatures were in good agreement with the corresponding experimental values.
Aspects that are ignored by most authors are the effect of convective heat flow due to material deformation and the asymmetry of heat generation due to the much higher pressure at the back of the shoulder. The former requires a prediction of the material flow around the tool, which is difficult to implement in most (nonfluid) solvers, which only predict weld temperature.
Recently, Assidi et al. in [11] used an ALE formulation implemented into the Forge3 software with a splitting approach and an adaptive remeshing scheme based on error estimation. In [21] the residual stresses in a 3D FE model were predicted for the FSW simulation of butt joints through a single block approach. The model was able to predict the residual stresses by considering only thermal actions. Buffa el al. [21] simulated the welding process using a continuous rigidviscoplastic finite element model in a single block approach through the Lagrangian implicit software, DEFORM3DTM. Then, the temperature histories extracted at each node of the model were transferred to another finite element model considering elastoplastic behavior of the material. The map of the residual stress was extrapolated from the numerical model along several directions by considering thermal actions only.
Santiago et al. [101] developed a simplified computational model taking into account the real geometry of the tool, i.e. the probe thread, and using an ALE formulation. They considered also a simplified friction model to take into account different slip/stick conditions at the pin shoulder/workpiece interface.
The more recent works performed in this direction are those presented in this paper. Agelet de Saracibar et al. ([3], [4]), Agelet de Saracibar et al. ([5], [6]), Chiumenti et al. [40], and Dialami et al. ( [50], [51], [52]) used a subgrid scale finite element stabilized mixed velocity/pressure/temperature formulation for coupled thermorigidplastic models, using Eulerian and Arbitrary Lagrangian Eulerian (ALE) formalisms, for the numerical simulation of FSW processes. They used ASGS and OSGS methods and quasistatic subgrid scales, neglecting the subgrid scale pressure and using the finite element component of the velocity in the convective term of the energy balance equation. Chiumenti et al. [40], Dialami et al. [50], and Chiumenti et al. [41] developed an apropos kinematic framework for the numerical simulation of FSW processes and compared with a solid approach in [22] and [24]. A combination of ALE, Eulerian and Lagrangian descriptions at different zones of the computational domain and an efficient coupling strategy was proposed. The resulting apropos kinematic setting efficiently permitted to treat arbitrary pin geometries and facilitates the application of boundary conditions. The formulation was implemented in an enhanced version of the finite element code COMET [30] developed by the authors at the International Centre for Numerical Methods in Engineering (CIMNE). Chiumenti et al. [41] used a novel stressaccurate FE technology for highly nonlinear analysis with incompressibility constraints typically found in the numerical simulation of FSW processes. They used a mixed linear piecewise interpolation for displacement, pressure and stress fields, respectively, resulting in an enhanced stress field approximation which enables for stress accurate results in nonlinear computational mechanics.
This paper presents the current state of the art in the numerical simulation of FSW processes. The outline is as follows:
In section 2 the thermal problem for the simulation of the FSW process at both local and global level is formulated. The thermal problem is governed by the enthalpy based balance of energy equation. Heat generation via viscous dissipation as well as frictional heating due to the contact is taken into account. Thermal convection and radiation boundary conditions are also considered.
In section 3 the mechanical problem for the simulation of the welding process is formulated. The mechanical problem is described by the balance of momentum equation.
The material behavior is either thermoelastoviscoplastic (global level analysis) or thermorigidviscoplastic (local level analysis).
Section 4 describes the discrete FE modeling and subgrid scale stabilization to deal with incompressibility and heat convection dominated problems. The multiscale stabilization method is introduced and an approximation of the subgrid scale variables together with the stabilization parameters is given. Algebraic Subgrid Scale (ASGS) and Orthogonal Subgrid Scale (OSGS) methods for mixed velocity, pressure and temperature linear elements are used. It is shown how the classical GLS and SUPG methods can be recovered as a particular case of the ASGS method.
Section 5 is devoted to the description of the proposed kinematic framework to simulate the FSW process. A novel numerical strategy to model FSW is presented. Using the ArbitraryLagrangianEulerian kinematic framework, the overall computational domain is divided into subdomains associating an apropos kinematic framework to each one of them. A combination of ALE, Lagrangian and Eulerian formulations for the different domain parts is proposed. Coupling between each domain is explained in detail including the friction contact. The strategy adopted to deal with an accurate definition of the boundary conditions is presented.
Section 6 compares a Solid Mechanics and a Fluid Mechanics approach for the numerical simulation of FSW processess. The main features of each one of those approaches including their main advantages and drawbacks are discussed.
Section 7 deals with the visualization of material flow during the welding process. A particle tracing technique is used, to visualize the trajectories of any material points. This method can be naturally applied to the simulation of the material flow in welding simulations. The trajectories of the material points are integrated from the velocity field obtained in the simulation at the postprocess level.
Section 8 is devoted to the description of the Localglobal strategy for the calculation of residual stresses in FSW process. The heat power is obtained at local level and then is transferred to the global one in order to compute the residual stresses.
Section 9 summarizes the work and provides a critical overview of the goals achieved in the simulation of FSW processes. Innovative features of the work are highlighted. Conclusions are drawn with regard to the fields of application and the intrinsic limits of the presented methods.
The governing equation representing the thermal problem is the balance of energy equation. This equation controls the temperature evolution and can be stated as ([1], [28], [38]):

(1) 
where , and are the density at the reference configuration, the specific heat, temperature, the volumetric heat source introduced into the system by plastic dissipation and the heat flux, respectively.
Depending on the aim of the analysis, whether it is local or global, there are different definitions for the heat source, :
Temperature field in a FSW process (Global level) 
Figure 6: Temperature field in a FSW process (Global level) 
Temperature evolution obtained from global level analysis compared with experiment 
Figure 7: Temperature evolution obtained from global level analysis compared with experiment 
Temperature contour field in FSW process (Local level). 
Figure 8: Temperature contour field in FSW process (Local level). 
Temperature evolution obtained from local level analysis compared with experiment 
Figure 9: Temperature evolution obtained from local level analysis compared with experiment 
At the global level, the Lagrangian framework is used to describe the thermal problem. According to this framework, the energy equation reads simply

(2) 
In global level analysis, the heat source is the (known) power input for the problem, obtained by solving the local level analysis.
Figure 6 shows the temperature field and moving heat source in a global level simulation.
The temperature field obtained from a global level analysis for different lines parallel to the weld line at the bottom surface is compared with those obtained from experiments is depicted in figure 7.
At the local level, where the focus of the simulation is the heat affected zone (fixed in the space together with the reference frame), the Eulerian framework is considered. According to this kinematic setting, the energy equation can be rewritten as

(3) 
where is the (spatial) velocity (see section 5 for further details).
In the FSW process, the heat source introduced into the system is due to the mechanical dissipation and it is computed as a function of the plastic strain rate, , and the deviatoric stresses, , as:

(4) 
where is the fraction of the total plastic energy converted into heat.
Application of the local level analysis in FSW process can be seen in Figure 8.
The temperature field obtained from a local level analysis for different lines parallel to the weld line at the bottom surface is compared with those obtained from experiments is depicted in figure 9.
Thermal boundary condition 
Figure 10: Thermal boundary condition 
Let us denote by an open and bounded domain. The boundary can be split into and such that , where fluxes (on ) and temperatures (on ) are prescribed for the heat transfer analysis as

(5) 
In figure 10, and represent the tool and the workpiece domains, respectively.
The initial condition for the transient thermal problem in terms of the initial temperature field is: .
On free surfaces, the heat flux is dissipated through the boundaries to the environment by heat convection, expressed by Newton's law as:

(6) 
where is the heat transfer coefficient by convection, is the surrounding environment temperature and is the temperature of the body surface.
Another heat dissipation mechanism is heat loss due to radiation. Heat radiation flux is computed using the StefanBoltzmann law:

(7) 
where is the Stefan–Boltzmann constant and is the emissivity factor.
Heat transfer by conduction involves transfer of energy within a material without material flow. The thermal constitutive model for conductive heat flux is defined according to the isotropic conduction law of Fourier. It is computed in terms of the temperature gradient, , and the (temperature dependent) thermal conductivity, , as:

(8) 
Velocity (top) and temperature fields (bottom) obtained with fully stick condition between pin and workpiece. 
Figure 11: Velocity (top) and temperature fields (bottom) obtained with fully stick condition between pin and workpiece. 
Adopting the classical friction model based on Coulomb's law, the so called slip function is defined as [7]:

(9) 
where is the friction coefficient, which can result in a nonlinear function of temperature and slip displacement . and are the normal and tangential components of the traction vector, , at the contact interface, respectively:

where is the unit vector normal to the contact interface.
This given, both stick and slip mechanisms can be recovered using the unified format:

(12) 
together with the KuhnTucker conditions defined in terms of the slip function, , and the slip multiplier, , as:

(13) 
where is a penalty parameter (regularization of the Heaviside function).
and are the tangential and the normal components of the total relative displacement, computed as:


(16) 
The normal traction vector is obtained with a further penalization as:

(17) 
where is the normal penalty parameter (not necessarily equal to ), which enforces nonpenetration in the normal direction.

(18) 
while the normal component of the traction vector is given by Eq. (17).
The relative velocity between two bodies in contact is the cause of heat generation by friction. This is one of the key mechanisms of generating heat in the FSW process. When the driving variable is the velocity field, , it is very convenient to use a Norton type friction law.
The tangential component of the traction vector at the contact interface, , is defined as:

(19) 
where is the equivalent friction coefficient. is the (temperature dependent) material consistency, is the strain rate sensitivity and is the tangential unit vector, defined in terms of the relative tangential velocity at the contact interface.
This given, the heat flux generated by Norton's friction law reads:

The total amount of heat generated by the friction dissipation is split into the fraction absorbed by the bodies in contact. The amount of heat absorbed by the first body, , and by the second body, , depends on the thermal diffusivity, , of the two materials in contact as:

The more diffusive is the material of one part in comparison with the other part, the more heat is absorbed by it.
Considering Eqs. (2) and (8), the weak form of the thermal problem at global level is defined over the integration domain and its boundary as:

where is the test function of the temperature field, while the thermal work, , is defined as:

(23) 
Note that the term due to friction () appears only at local level.
Similarly, considering Eqs. (3) and (8), the weak form of the problem at local level is written as

In the framework of the standard Galerkin finite element method, the discrete counterpart of the weak form for the thermal problem can be written as


where is the discrete temperature field.
Mechanical results in a FSWprocess (local analysis). 
Figure 12: Mechanical results in a FSWprocess (local analysis). 
In both cases, assuming quasistatic conditions, the local form of the balance of momentum equation, also known as Cauchy's equation of motion, is given by

(27) 
where is the Cauchy's stress tensor, is the body forces vector per unit of mass, is the density in the reference configuration and is the divergence operator.
In the case of metallic materials where the behavior of inelastic (plastic) strains is isochoric (plastic deformations are deviatoric), it is convenient to split the stress tensor, , into volumetric and deviatoric parts:

(28) 
where is the pressure and is the stress deviator. Similarly, the strain tensor, can be split into volumetric, and the deviatoric, parts.

(29) 
Therefore, the local form of the mechanical problem can be stated as:

(30) 
where and are defined according to the constitutive equations selected for the local or global analysis, respectively.
Velocity streamlines around a triangular pin (left) and a trivex pin (right) 
Figure 13: Velocity streamlines around a triangular pin (left) and a trivex pin (right) 
Plastic dissipation under a FSWtool. 
Figure 14: Plastic dissipation under a FSWtool. 
Longitudinal stress (XX) field in a global level simulation of FSW after cooling. 
Figure 15: Longitudinal stress (XX) field in a global level simulation of FSW after cooling. 
Figure 12 shows the results of the mechanical simulation at the local level. Pressure, velocity, strainrate and stress fields are presented around the rotating triangular pin. The peaks in mechanical variables around the pin can be seen. Due to the high rotational velocity of the pin, the asymmetry of the results are not visible.
The streamlines and velocity vectors around the triangular pin are shown in figure 13. One can see the strong effect of rotational velocity and the differences between retreating and advancing sides.
Plastic dissipation under the FSW tool with cylindrical pin is illustrated in figure 14. Note that plastic dissipation is one of the main sources of heat in FSW. Effect of the shoulder is clearly important and cannot be neglected. This obliges modelers to perform simulation in 3D in order to obtain accurate results.
Figure 15 shows the results of the mechanical simulation at global level (residual stresses). The figure illustrates the effect of a moving heat source on the stress XX field after fixture release and cooling phase.
According to the general perspective of the paper, which studies the FSW process at both local and global level, this section is divided into two parts: the first part describes the mechanical constitutive equations used at global level; the second part is devoted to the description of the mechanical constitutive model used at the local level.
When the effects induced by the process on the structural behavior are the target of the study, the global level analysis is preferred. These effects are distortions, residual stresses or weaknesses along the welding line, among others and the displacement field is the driving variable.
At global level, the stress deviator, , is generally defined through the constitutive equation as a function of displacements and temperature The heat source moves along the welding path and the temperature range varies from room temperature to very high temperatures close to the melting point. Consequently, the material behavior changes from elastoplastic (at room temperature) to pure viscous (close to the melting point). This evolution can be defined by a thermoelastoviscoplastic constitutive model. Close to the melting point, the viscous behavior dominates and the elastic limit gradually decreases.
To this end, it can be written:

where is the bulk modulus, which controls the material compressibility and is the shear modulus. Deformations and are the volumetric thermal deformation, the elastic volumetric strain, the viscoplastic and the elastic deviatoric strain, respectively. The thermal deformation is expressed as

(33) 
where is the thermal expansion coefficient.
Temperature dependent material properties of 304L stainless steel. 
Figure 16: Temperature dependent material properties of 304L stainless steel. 
Figure 16 shows the temperaturedependent material property for 304L stainless steel. The mechanical material parameters can be obtained by fitting the constitutive model with uniaxial, tensile stressstrain curves obtained at different temperatures.
The viscoplastic deformation is defined through appropriate evolution laws

where is the viscoplastic yield surface, is the unit normal to the yield surface and is the viscoplastic parameter, when a ratedependent evolution law is assumed. is the plastic viscosity associated to the viscoplastic model.
As in this model the elastic effects are not neglected, the residual stresses can be evaluated.
Note that in the liquidlike phase, the viscoplastic potential

(36) 
degenerates into a purelyviscous potential , when the coefficient . As a consequence, considering the definition of the viscoplastic multiplier, it transforms into . This means that the viscoplastic model reduces to a NortonHoff model which is a viscous relationship between the deviator of the stress tensor, , and the viscoplastic strain rate, , governed by the viscous parameter, .

(37) 
When the heat affected zone is the focus of the study, the local level analysis in Eulerian format is preferred. The stress deviator is defined as a function of temperature, and the velocities, rather than displacements. The constitutive laws are typically formulated in terms of strainrate rather than strain. The FSW process is characterized by very high strain rates as well as a wide temperature range going from the environmental temperature to the one close to the melting point. Hence, instead of a thermoelastoviscoplastic model (generally adopted for metals in the Lagrangian formulation), a thermorigidviscoplastic behavior is usually introduced within an Eulerian/ALE framework. The elastic strains are neglected compared to the viscoplastic flow and, therefore, precluding the computation of residual stresses.
When the problem is defined in terms of velocities, the strain rate tensor is obtained as . The split into volumetric and deviatoric parts results in:

(38) 
where and are the volumetric and the deviatoric parts of the strain rate tensor, respectively.
The incompressibility constraint is often adopted by assuming that the volumetric changes, including thermal deformation, are negligible in comparison with elastoplastic (deviatoric) deformation:

(39) 
so that .
According to the split of the stress tensor introduced in Eq. (28), different ratedependent constitutive models can be used for modeling the welding process. The constitutive models are usually defined by the relationship between the deviatoric parts of the stress and the strain rate, . This reduces into the definition of viscoplastic constitutive models of the form:

(40) 
where is the effective viscosity coefficient of the material. The following additive decomposition of the deviatoric strain rate is assumed:

(41) 
where and are the elastic and viscoplastic parts, respectively. In local level FSW, the elastic part of the strain tensor, , is negligible compared with the viscoplastic component, , so that:

(42) 
This means that the total deformation is viscoplastic and a rigidviscoplastic model is recovered. The rigidviscoplastic models can be interpreted as nonNewtonian laws where the viscosity is a nonlinear function of the strainrate Typical constitutive equations used for nonNewtonian modeling are listed in Table constitutive lists. The coefficients and in these models represent the asymptotic values of the viscosity at zero and infinite strainrates, respectively. In the Newtonian case
In the present paper, depending on the definition for , three models are introduced from the groups (b), (c) and (f): NortonHoff, SheppardWright and Carreau constitutive models.
NortonHoff model The NortonHoff model was originally introduced by Norton [90] in order to describe the unidimensional creep of steel at high temperature, and extended by Hoff [70] to multidimensional situations. This model can be easily characterized by simple mechanical tests and accurately describes the behavior of materials with significant viscous and creep response to loading. The NortonHoff model is best suited for the materials that are tough and highly rate sensitive. According to this constitutive model, the effective viscosity is expressed as a function of the temperature and the equivalent plastic strain rate, , as:

(43) 
where and are the (temperature dependent) rate sensitivity and viscosity parameters, respectively (Figure 17). Note that:

(44) 
where defines the plasticflow direction.

(45) 
that is, a linear relationship between stresses and strainrates where the constant of proportionality is the viscosity parameter.
The great advantage of the NortonHoff model is its simplicity as it only depends on two (temperature dependent) parameters. For the majority of metals, the ratesensitivity parameter is usually within the range which results in a very nonlinear (nonNewtonian) material behavior (e.g. aluminium and steel alloys). The nonlinearity of the NortonHoff model increases as diminishes. To deal with this nonlinearity, can be computed in a series of steps starting from the value in the first step and reaching its real value at step . This implies that at each timestep the problem is solved with constitutive rate sensitivity :

(46) 
Stress vs. strain rate affected by rate sensitivity variation. 
Figure 17: Stress vs. strain rate affected by rate sensitivity variation. 
SheppardWright model The equation used to describe the flow strength of alloys was introduced initially by Sellars and Tegart [105]. Sheppard and Wright [106] introduced an alternative rigid viscoplastic formulation extending this formulation. According to the SheppardWright model, the effective viscosity coefficient is defined in terms of effective flow stress (i.e. norm of deviatoric stress) and effective strain rate (i.e. norm of deviatoric strain rate) as

(47) 
where the effective stress depends on the strain rate and the temperature field (below the metals solidus temperature) and it is defined as

(48) 
Effect of activation energy on stress vs. strain rate. 
Figure 18: Effect of activation energy on stress vs. strain rate. 
where , the Zener–Hollomon parameter, represents the temperature compensated effective strain rate [121]

(49) 
is the universal gas constant and is the absolute temperature (in Kelvin). The material constants and the apparent activation energy, are derived by fitting experimental response of the material, and are all independent of temperature (Figure 18). The advantage of the SheppardWright model in comparison with the NortonHoff model is the possibility of a better calibration of the material behavior in the entire temperature range from the environment temperature to the melting point.
Carreau model According to Carreau's constitutive model ([26], [102]), the effective viscosity coefficient is defined as

(50) 
StressStrain rate behavior of a Carreau model with temperature change. 
Figure 19: StressStrain rate behavior of a Carreau model with temperature change. 
where , , and are the saturation viscosities at zero and infinite limits, the plastic strain rate corresponding to the yield limit and the exponent, respectively. The yield limit is defined by a Johnson–Cook power of the form

(51) 
where , , and are melting temperature, room temperature, yield limit at room temperature and Johnson–Cook exponent, respectively. The main advantage of this constitutive model is that the range of admissible viscosity values is bounded by the viscosity at room temperature and for the liquid phase. These limit bounds will resolve the problem related to singularity of the nonNewtonian constitutive models and thus use of numerical tricks such as definition of cutoff values on the strain rates. Figure 19 shows the stressstrain rate behavior of a Carreau model for different temperatures. Neither Carreau and SheppardWright nor Norton Hoff model consider thermoelastic strains. Therefore the proposed constitutive laws cannot predict residual stresses after joining and cooling back to the room temperature. The prediction of residual stresses is the outcome of the process study at global level where the whole structure is simulated assuming moving heat source obtained from the local based simulation.
The mixed formulation is the most classical formulation to tackle the incompressible limit, where both displacements and pressures are used as master fields. It is often used for modeling nearly incompressible (Poisson's ratio above 0.49) or completely incompressible material behavior (Poisson's ratio equal to 0.50). The standard irreducible formulation is simply not suitable for incompressible conditions. The mixed formulation solves the balance of momentum equation (30) and the volumetric part of the constitutive Eq. ( 31) in the weak form. Taking into account that is defined as a function of displacement, the weak form of the mixed mechanical problem at global level is defined over the integration domain and its boundary as:

(52) 
where and are the test functions of the displacement and pressure fields, respectively, while the mechanical work is defined as:

(53) 
Definition of the problem in the mixed format is not obligatory, however it is convenient for including both the compressible and incompressible cases. If and the bulk modulus such as in the liquidlike phase, the incompressibility condition is recovered: . In this case, the second Eq. (52) enforces (in weak sense) the incompressibility condition as a purely kinematic constraint. Similarly, according to Eqs. (30) and (39), and taking into account that is defined as a function of the velocity field, the weak form of the mixed mechanical problem at local level is defined over the integration domain and its boundary as:

(54) 
where is the test function of the velocity field.
For the global level analysis, the boundary is split into and , being such that the tractions are prescribed on while displacements are specified on . The corresponding Neumann and Dirichlet boundary conditions on and are defined as

(55) 
where and are the prescribed traction and velocity field, respectively and is the unit outward normal to the surface . Similarly, for the local level analysis, the boundary is split into and , being such that the velocities are specified on and the corresponding Dirichlet and Neumann boundary conditions on and are defined.
In the framework of the standard Galerkin finite element method, the discrete counterparts of the weak forms for the mechanical problem can be written as

(56) 

(57) 
where and are the discrete counterparts of the respective fields.
In this section, a stabilized finite element formulation based on the Variational Multiscale (VMS) stabilization method is adopted. This stabilization method allows the use of equal order interpolation for the velocity and pressure fields, bypassing the need to satisfy the infsup condition (see first subsection), and avoiding the oscillations arising in convective dominated problems (see second subsection).
It is known that the incompressibility constraint or under isochoric strain conditions, the use of an irreducible velocitybased (or displacementbased) formulation leads to pressure locking situations. It must be pointed out that when using a J2plasticity model (in solid mechanics) or one of the nonNewtonian laws presented in the previous section (for the Eulerian framework), it is usual to get an isochoric strain field, that is, most of the deformations are deviatoric and the volumetric part is negligible. The use of a mixed formulation is a wellknown remedy. However, when the mixed Galerkin formulation is used, compatible spaces for the pressure and the velocity field have to be considered and they must satisfy the LadyzhenskayaBabuska–Brezzi (LBB) stability condition [18]. The LBB compatibility condition restricts the choice of finite element spaces and for velocities (or displacements) and pressure, respectively. It states that in order to obtain a stable formulation, and must be chosen such that the following infsup condition is satisfied:

(58) 
where parameter is independent of the mesh size . If the LBB compatibility condition is satisfied, there exists a unique solution and a (determined up to an arbitrary constant in the case of purely Dirichlet boundary conditions). Unfortunately, standard Galerkin mixed elements with continuous equal order linear velocity/pressure interpolation violate the LBB condition, leading to instabilities in the pressure field and poor numerical performance. Stable formulations can be obtained either by choosing velocity/pressure pairs that satisfies the LBB condition (e.g. P2/P0 [71]) or by using pressure stabilization methods. Using LBB stable elements, such as the Q2/Q1 (TaylorHood element, based upon continuous quadratic velocity, continuous bilinear pressure) or the mini element (continuous linear velocity+bubble function, continuous linear pressure) is difficult in practice, due to the complexity of the associated implementation and lower computational efficiency than that of simplicial elements. Fortunately, there exist possibilities for circumventing the LBB condition, which permit the use of equal order velocity/pressure interpolations, by modifying the original Galerkin problem, adding stabilization terms. Essentially, in the majority of stabilization methods, the stabilizing terms added to the original Galerkin formulation involve the residual of the balance of momentum equation. These methods are called residualbased. A stabilization method is consistent if, on convergence, the solution of the modified problem is exactly the solution of the original problem. The subgrid method was originally introduced in [73] as a powerful technology known for its enhanced accuracy properties. The exact solution is split into two parts, corresponding to different scales or levels of resolution: the one that can be captured by the finite element mesh (coarse) and the other one, called the subgrid (finer) part. While the first one can be resolved by the finite element method, the latter one cannot. The particular approximation used for the subgrid scale defines the numerical strategy. The solution of the continuous problem contains components from both scales. It is necessary to, somehow, include the effect of both scales in the approximation to get the stable solution of the discrete problem. The aim is to find an approximate solution for the subgrid scale and to include it into the discrete finite element solution. The effect of this finer scale can be included, at least locally, to enhance the stability of the pressure in the mixed formulation ([Chiumenti et al. 2002], [2], [37], [Cervera et al. 2003]). To this end, the solution space, , is split into two parts:

(59) 
denotes the finite element solution while , is the one that completes in and it is called the subgrid scale ([73] and [72]). Accordingly, the velocity and pressure fields of the mixed problem are expressed as (Figure 20)

Finite element and subgrid scale parts of the solution U=vpmath 
Figure 20: Finite element and subgrid scale parts of the solution U=vpmath 
where is the velocity field of the finite element scale and is the enhancement of the velocity field belonging to the subgrid scale. Note that no subgrid scale contribution is considered for the pressure field. It is reasonable to assume that the subgrid velocities are sufficiently small compared to . They can be viewed as a high frequency perturbation of the finite element field, which cannot be resolved in . It can also be assumed that and its variation vanish on the boundary . The stress deviator can be similarly decomposed into two different parts:

(61) 
Notice that the stress is not a linear function of and, therefore, taking and is a further approximation. Different choices are available for the approximation of the velocity subgrid scale. Two possible options correspond to the ASGS and the OSGS methods. Using ASGS method, is approximated as , where is the residual of the momentum equation in finite element space.

(62) 
where vanishes when using linear triangular or tetrahedral elements and the body forces, , are negligible. This given:

(63) 
After substituting into the weak form of both the balance of momentum equation and the constitutive equation, the discrete ASGS stabilized weak form of the mechanical problem for linear elements simplifies to :

(64) 
Representation of OSGS solution space 
Figure 21: Representation of OSGS solution space 
The orthogonal projection of a variable can be computed as: . The resulting discrete OSGS stabilized weak form of the mechanical problem is:

(65) 
where and are the smooth projection of the pressure gradient and its variation onto the finite element space, respectively. Observe that the stabilization term introduced in the second equation of (65) is computed as a function of the orthogonal projection of the residual of the momentum balance equation as:

Pressure field obtained using ASGS (top) and OSGS (bottom) stabilization method. 
Figure 22: Pressure field obtained using ASGS (top) and OSGS (bottom) stabilization method. 
Streamlines obtained using ASGS (top) and OSGS (bottom) stabilization method. 
Figure 23: Streamlines obtained using ASGS (top) and OSGS (bottom) stabilization method. 
The stabilization parameter , is computed as:

(66) 
where is the element size, is a constant and is the parameter defined according to the constitutive model (either effective viscosity or shear modulus). is usually temperature dependent and therefore the stabilization term is a temperature dependent parameter. Figure 22 compares results of ASGS and OSGS pressure stabilization methods in a parallel flow passing through a tube. The results obtained using OSGS are nodally exact, and totally unaffected by the boundary conditions [6]. Figure 23 shows the streamlines. The results obtained using OSGS are totally unaffected by the boundary conditions.
When using an ALE or Eulerian framework for the description of the problem, a convective term coming from the material time derivative arises in the equations. It is wellknown that in diffusion dominated problem, the solution of the balance of energy equation is stable. However, in convection dominated problem, the stabilizing effect of the diffusion term is insufficient and instabilities appear in the solution. The threshold between stable and unstable conditions is usually expressed in terms of the Peclet number. The typical Peclet number for FSW processes ranges from to . For this range of the Peclet number, the effect of the convective term cannot be neglected and the solution of the thermal problem is found to be instable. According to the subgrid scale method, the solution space of the temperature field is split into two parts: the finite element space and the space for the subgrid scale.

(67) 
Therefore, the temperature field can be approximated as

(68) 
where is the temperature component of the (coarse) finite element scale and is the enhancement of the temperature field corresponding to the (finer) subgrid scale. Let us also consider the corresponding variations and , respectively. Also in this case, we consider two possible choices for the approximation of the temperature subgrid scales corresponding to the ASGS and the OSGS methods. In the ASGS method, is approximated as , where is the residual of the energy balance equation in the finite element scale.

(69) 
where the diffusive term, vanishes when using linear triangular or tetrahedral meshes. In this case, the discrete thermal equation for linear elements can be stabilized by adding a residual based ASGS stabilization term of the form:

(70) 
to the Galerkin weak form of the problem. This given:

(71) 
It can be observed that for linear elements, the ASGS method coincides with the SUPG method. As an alternative, the space orthogonal to the finite element space can be adopted for the approximation of the unknown subgrid space . This means that the solution space is approximated as . According to the OSGS method, is defined as . Taking into account that , the resulting discrete OSGS stabilized weak form of the thermal problem for linear elements simplifies to:

(72) 
where is the smooth projection of the convective term.

(73) 
The mesh dependent stabilization parameter defined for the thermal (convective) problem is

(74) 
Comparison between the temperature distribution at different times (10, 40 and 70) and at the steady state. 
Figure 24: Comparison between the temperature distribution at different times (10, 40 and 70) and at the steady state. 
Comparison between the temperature distribution at the center line at different times (10, 40 and 70) and at the steady state, obtained using GLS pressure stabilization method and OSGS and SUPG convection stabilization methods are shown in figure 24 in a parallel flow passing through a tube.
The kinematic framework used for the simulation of FSW process is different if a local or global analysis is performed. When the process is simulated at global level, the Lagrangian framework is preferred. In this case the movement of the material points, , coincides with the mesh nodes. The material and local time derivative of the temperature concide:

(75) 
If the process is studied at the local level, ALE or Eulerian frameworks are a convenient choice. In this case, the material derivative depends on the local derivative plus a convective term of the form:

(76) 
where is the spatial velocity and is the spatial gradient. The convective term accounts for the movement of the material points with respect to a fixed mesh. In the ALE framework, the reference system neither moves together with the material (Lagrangian) nor is fixed (Eulerian) but it can move arbitrary. In this case, the material derivative depends on the local derivative plus a convective term as a function of the relative velocity . It is defined as

(77) 
where the relative velocity is . and are the material and the mesh velocities, respectively . In this case, the movement of the mesh must be calculated at each timestep. In general, an adhoc methodology is necessary to compute the position of the mesh at each timestep of the analysis. This is known to be one of the main difficulties of the ALE method. At this point, it is convenient to introduce a feasible strategy and an apropos kinematic framework for the simulation of the FSW process. In the FSW process, the pin rotates at constant angular velocity and, at the same time, advances with a constant speed. For the sake of convenience, in the numerical simulation, the movement of the pin is defined by the pure rotation around its axis while the advancing speed (in the opposite direction) is imposed to both the workpiece and the backplate. Taking into account that, during the welding process, the pin is rotating with a very high speed (e.g. 505000 RPM, depending on the workpiece material), a fully Lagrangian approach (which follows the material particles of the continuum in their motion) leads to a computational overhead. The reason is that in a Lagrangian analysis, the mesh is attached to the material and thus the mesh deforms with it. The material in the stirring zone suffers very large deformations at high strainrates. As a consequence, a fine mesh and a continuous remeshing are required to avoid an excessive mesh distortion. This usually leads to computationally intensive analyses, as well as, to a general loss of solution accuracy due to the interpolation process necessary to move both nodal and Gaussian variables from mesh to mesh. Application of a fully Eulerian approach is feasible exclusively for pins with a circular crosssection. When the pin is not cylindrical, the boundaries of the model are continuously changing according to the actual position/rotation of the pin. As a consequence, the integration domain must be redefined at each timestep of the simulation. In modeling FSW, the choice of kinematic frameworks for different domain zones (regions) crucially impacts on the computational efficiency and the solution quality. For defining the suitable frameworks several observations must be made:
The geometry discretization 
Figure 25: The geometry discretization 
In the proposed methodology, the original geometry is divided into three distinct zones (Figure 25) which are: the pin, the workpiece and the stirzone (process zone close to the pin where most of the material deformation takes place). An apropos kinematic framework is adopted for each one of these zones.

(78) 
where is the angular velocity of the pin and is the position of any grid point with respect to the rotation axis.
The ALE kinematic framework can recover the Lagrangian and the Eulerian frameworks as limit cases. Denoting by the mesh velocity, the convective term can be written as , the convective term vanishes in the Lagrangian framework ( ) and simplifies to in the Eulerian zone ( ). Table 2 summarizes the computational framework together with the solution hypotheses for the pin, the workpiece and the stirzone. According to the geometry subdivision and kinematic frameworks introduced, two different types of domain interactions exist. The coupling between each part needs a special attention: on one hand, the connection between the fixed mesh of the workpiece and the ALE mesh used for the stirzone must be considered. On the other hand, the frictional contact between the stirzone and the pin (see section 2.5) must be defined. The coupling at the interface between the ALE (stirzone) and the Eulerian (workpiece) parts requires a special treatment as the mesh is fixed on one side (Eulerian) and it moves on the other side (ALE). Physically, workpiece and stirzone are not separated, they belong to the same metal sheet even if in the numerical model, a relative movement of the two computational subdomains exists. Therefore, the objective is to get continuous fields for all the state variables and at the interface between the two zones. The coupling is performed using a nodetonode link or nodetoface contact approach.Nodetonode link approach. 
Figure 26: Nodetonode link approach. 
In the first and simpler case, each node of the mesh on the interface is duplicated (Figure 26). One node belongs to the ALE region and the other one to the Eulerian region. At every mesh movement, for a given node of the ALE surface, the corresponding node of the Eulerian side is found and a link between the two nodes is created. A search algorithm restricted to the interface nodes is performed in order to identify all the nodetonode connections at the interface. This leads to an adhoc assembling procedure, where the contributions (elemental residuals and tangent matrices) of the stirzone nodes are summed to the corresponding contributions of the workpiece nodes. This approach requires a coincident mesh at the interface. The timestep is chosen such that the two meshes (ALE and Eulerian) are synchronized to keep the Courant number, at the sliding interface. The nodetonode linking approach is simple to achieve for 2D analysis but it becomes complex when an automatic 3D mesh generator (which usually supports unstructured tetrahedral meshes) is used. This approach can be applied in 3D analysis if a structured mesh is defined at the interface surfaces. It should be meshed in such way that at each timestep, the ALE mesh at the interface slides from one node to another. This requires coincident and equispaced mesh at the contact surfaces. If the surface meshes are not coincident, a nodetoface approach can be used. For this more complex approach, the mesh is not sychronized. Therefore, at the sliding interface, contact elements are created. At the contact interface, each node of the stirzone is projected on the workpiece surface. Once the contact elements are generated, whether the Lagrange multiplier method, the augmented Lagrange method, or the penalty method can be used to apply the constraint between the problem variables, , , . In the penalty method, very large values of both stiffness and thermal resistivity representing the penalty parameters are assigned to the contact elements. This is performed to ensure the most rigid/conductive link between the two domains in both the mechanical and the thermal problems. In this case, the integration domains are connected by forcing coincident heat fluxes and normal tractions on both sides of the contact surface. In this approach, at each timestep, the contact elements need to be regenerated according to a (rather timeconsuming) closestpointprojection algorithm.
As it was pointed out in Section 5, the most suitable kinematic framework for the numerical simulation of FSW processes may be different if local or global level analysis is performed. When the numerical simulation of FSW processes is performed at global (structural component) level, considering the full component to be joined, the use of a Lagrangian kinematic description, within a Solid Mechanics approach, is the most suitable option. On the other hand, when the numerical simulation of FSW processes is performed at local level, the use of an adhoc Apropos Kinematic Framework (AKF) which makes use of an efficient combination of Lagrangian (tool), Eulerian (workpieces) and ALE (stirring zone) descriptions for the different computational subdomains, within either a Solid Mechanics or a Fluid Mechanics approach, is the most suitable option ([22], [23], [24], [50]). In this section a comparison between a Solid Mechanics and a Fluid Mechanics approach for the numerical simulation of FSW processess at local level will be performed ([22], [23], [24]). The main features of each one of those approaches, pointing out which are their main advantages and drawbacks, will be presented.
In the Solid Mechanics approach, the governing equations of the coupled thermomechanical problem are written in terms of the displacement vector field and temperature scalar field, which are taken as primary variables ( [22], [23], [24]). The original coupled thermomechanical problem is solved using a staggered procedure which arises from an isothermal operator split. The isothermal operator split defines a mechanical problem (at constant temperature) and a thermal problem (at fixed mechanical configuration, keeping frozen the current configuration). The primary variable of the mechanical problem (mechanical variable) is the displacement vector field, while the primary variable of the thermal problem (thermal variable) is the temperature field. The isothermal operator split applied to the governing equations of the coupled thermomechanical problem leads to a Product Formula Algorithm (PFA), where the original coupled thermomechanical problem is solved by means of a staggered timemarching scheme, where the mechanical and thermal subproblems are solved sequentially. At each time step, an isothermal mechanical subproblem is solved first, keeping constant the temperature computed at the previous time step. Then, the thermal subproblem is solved next, keeping constant the configuration computed in the mechanical subproblem.
In the Fluid Mechanics approach, the governing equations of the (fully incompressible) coupled thermomechanical problem are written in mixed form in terms of the velocity vector field, pressure scalar field, and temperature scalar field, which are taken as primary variables ([22], [23], [24], [50]). The original coupled thermomechanical problem is solved using a staggered procedure which arises from an isothermal operator split. The isothermal operator split defines a fully incompressible mechanical problem (at constant temperature) and a thermal problem (at a fixed mechanical configuration, keeping frozen the current velocity and pressure fields). The primary variables of the mechanical problem (mechanical variables) are the velocity and the pressure fields, while the primary variable for the thermal problem (thermal variable) is the temperature. The isothermal operator split applied to the governing equations of the coupled thermomechanical problem leads to a Product Formula Algorithm (PFA), where the original coupled thermomechanical problem is solved by means of a staggered timemarching scheme, where the mechanical and thermal subproblems are solved sequentially. At each time step, a thermal subproblem is solved first, keeping constant the mechanical variables, velocities and pressures, computed at the previous time step. Then, an isothermal mechanical subproblem is solved next, keeping frozen the temperature field computed in the thermal subproblem.
In the Solid Mechanics approach, a standard Galerkin linear displacement/temperature finite element formulation is used. Quadrilateral (2D) or hexahedral (3D) finite elements are used for the domain discretization. The position and temperature fields are computed at each node of the elements. The stress tensor field and the internal variables are computed at each integration point of the elements. Four (2D) or eight (3D) integration points are used in each element. To overcome the locking phenomenon linked to the quasiincompressibility constraint, the pressure is considered constant over the element and it is computed only at a central integration point ([22], [23], [24]).
In the Fluid Mechanics approach, a subgrid scale stabilized mixed linear velocity/pressure/ temperature finite element formulation is used. The Orthogonal Subgrid Scale (OSS) stabilization method is used to solve both the pressure instability arising from the incompressibility constraint, and the instabilities arising from the convective term in the energy balance equation in advection dominant problems. Triangular (2D) or tetrahedral/hexahedral (3D) finite elements are used for the domain discretization. The velocity, pressure and temperature fields are computed at each node of the elements. The deviatoric component of the stress tensor field and the internal variables are computed at each integration point of the elements ([22], [23], [24], [50], [6]).
Computational subdomains of the plates: Stirzone (zone 1), transition zone (zone 2), and workpieces zone (zone 3) 
Figure 27: Computational subdomains of the plates: Stirzone (zone 1), transition zone (zone 2), and workpieces zone (zone 3) 
In the furthest zone from the tool (zone 3), a Eulerian description is a suitable option. Thus, in this region the mesh is fixed. This can be considered as a particular case of an ALE description where the mesh velocity is set to zero. The circular (2D) or cylindrical (3D) crown region connecting zones 1 and 3 can be viewed as a transition zone (zone 2). Two different techniques are proposed to deal with the numerical simulation of the FSW process in the transition zone. The solution scheme of the governing equations using an ALE description in the zones 1 and 3, as well as the treatment of the transition zone 2, is different in the Solid Mechanics and Fluid Mechanics approach. The tool is modeled as a rigidthermal body and a Lagrangian description is used in both the Solid Mechanics and Fluid Mechanics approaches.
In the Solid Mechanics approach, an ALE description is applied to each one of the three computational subdomains defined in the Figure 27. In the stirzone 1, the mesh velocity is prescribed to the rotational velocity of the tool. In the workpieces zone 3, the mesh velocity is prescribed to zero. Thus, in the zone 3, the mesh is fixed in the space, resulting in an Eulerian description, which is treated as a particular case of the ALE formulation. In the transition zone 2, the mesh velocity is prescribed using a linear interpolation of the prescribed mesh velocities on the boundaries of the computational subdomain, the rotational tool velocity at the interface with the zone 1, and zero velocity at the interface with the zone 3 ([22], [23], [24]). Note that the quality of the mesh in the zones 1 and 3 does not changes during the simulation, while the quality of the mesh in the zone 2 will quickly deteriorate. As the mesh distortion grows with time, a remeshing operation is periodically required. For one full rotation of the tool, the remeshing process is performed 30 times. The remeshing operation can be divided into two steps. First a better suited mesh is created. The relatively simple geometry of the transition zone allows an easy generation of a new quadrangular mesh (3D). Then, to carry out the computation over the new mesh, the state variables from the old mesh (the mesh used before the remeshing operation was applied) have to be transferred to the new one. Each field used to define the equilibrium state is transferred independently from the other ones. The data transfer method used is the Finite Volume Transfer Method (FVTM) with linear reconstruction of the fields [25]. Using an ALE description, the governing equations of the coupled thermomechanical problem, given by the mass continuity equation, the momentum balance equation and the energy balance equation, take the form:

where is the spatial mass density, is the specific heat capacity, is the material velocity vector field, is the convective velocity vector field, defined as the difference between the material velocity and the mesh velocity vector fields, is the Cauchy stress tensor field, is the body forces vector field per unit of mass, is the temperature field, is the mechanical dissipation rate per unit of spatial volumen, is the heat flux vector per unit of spatial surface, and is the mesh coordinate system of reference for the ALE description. To simplify the solution procedure and remain competitive against Lagrangian based formulations, the ALE system of equations (79) is solved using an operator split procedure, which yields to a product formula algorithm (PFA) ([22], [23], [24], [16], [17], [25], [53]). First, for each time step, a Lagrangian step is solved. This Lagrangian step is formally obtained setting the convective velocity to zero. During this Lagrangian step the mesh sticks to the material and, thus, the mesh velocity is equal to the material velocity. An iterative procedure is performed until a Lagrangian equilibrium configuration is obtained. The spatial weak form of the momentum and energy balance equations to be solved in this Lagrangian step takes the form:

where and are the test functions for the displacement vector field and temperature field, respectively, is the prescribed traction vector field on the boundary , is the prescribed normal heat flux on , is the conduction heat fluxes on , and are the convection and radiation normal heat fluxes on , respectively, and is the isotropic thermal conductivity parameter. The upper dots stand for material time derivatives. After the solution of the Lagrangian step has been obtained, a second step, denoted as Eulerian step, is performed. This Eulerian step is divided into two substeps. First the nodes of the mesh are relocated to a more suitable position given in terms of the mesh velocity. Thus, a new mesh with the same topology is obtained. Then, the nodal and internal variables at the integration points are mapped from the old mesh to the new one.
In the Fluid Mechanics approach, the transition zone 2 is collapsed to a zero thickness circular (2D) or cylindrical (3D) crown interface. Each node of the mesh of this interface is duplicated. One node is linked to the stirzone 1, where an ALE formulation is used, and the other one is linked to the workpieces zone 3, where a Eulerian formulation is used. The coupling between the two zones, zone 1 and zone 3, can be ideally performed using a specific nodetonode (NTN) direct contact link approach. This approach requires a coincident and equispaced mesh at the interface. At every mesh movement step, for each node of the zone 1 of the common interface, the corresponding node of the zone 3 of the common interface is found and a NTN link is created. Then, the boundary conditions and properties of the nodes of the zone 3 are copied to the corresponding linked nodes of the zone 1. The time step can be conveniently chosen such that the two interface meshes (ALE in zon 1 and Eulerian in zone 3) are always synchronized to keep the Courant number . In this case the ALE mesh would slide precisely from one Eulerian interface node to the next one at each time step. If the meshes at the interface are not coincident, or if they are coincident but not equispaced, the meshes cannot be synchronized and a more general nodetosegment (NTS) contact approach should be used. Further details can be found in Section 5. In the Fluid Mechanics approach an Apropos Kinematic Framework (AKF) which makes use of an efficient combination of Lagrangian (tool), Eulerian (workpieces, zone 3) and ALE (stirring zone, stir zone 1) descriptions for the different computational subdomains, is used. The ALE formulation used in the Fluid Mechanics approach is not based on an operator split, as it was the case for the ALE formulation used in the Solid Mechanics approach. The system of ALE equations, for each time step, is solved here directly, in a single step, including the convective terms written in terms of the convective velocity. Remarkably, remeshing and data mapping operations needed in the Solid Mechanics approach, are fully avoided in the Fluid Mechanics approach. The following hypothesis are introduced within the Fluid Mechanics approach ( [22], [23], [24], [50], [6]):
Taking into account the above hypothesis, using an ALE description and a mixed velocity/pressure formulation to avoid the locking due to the incompressibility constraint, the strong form of the governing equations in the Fluid Mechanics approach takes the form:

where is the reference mass density, is the specific heat capacity, is the material velocity vector field, is the convective velocity vector field, defined as the difference between the material velocity and the mesh velocity vector fields, is the deviatoric part of the Cauchy stress tensor field, is the pressure field, is the body forces vector field per unit of mass, is the temperature field, is the mechanical dissipation rate per unit of spatial volumen, is the heat flux vector per unit of spatial surface, and is the mesh coordinate system of reference for the ALE description. The weak form of the mass continuity for an incompressible material, momentum for a quasistatic problem, and energy balance equations takes the form:

where

Evolution of the temperature computed by the Solid Mechanics and Fluid Mechanics approaches at three different points. (a) Different zones of the model with the initial position of three selected points; (b) Temperature at point 1 (rotating with the tool); (c) Temperature at point 2 (rotating with the tool); (d) Temperature at point 3 (fixed in space) Bussetta2013 
Figure 28: Evolution of the temperature computed by the Solid Mechanics and Fluid Mechanics approaches at three different points. (a) Different zones of the model with the initial position of three selected points; (b) Temperature at point 1 (rotating with the tool); (c) Temperature at point 2 (rotating with the tool); (d) Temperature at point 3 (fixed in space) Bussetta2013 
In both the Solid Mechanics and Fluid Mechanics approaches, the tool is modeled as a rigid body, using a thermorigid constitutive model, and the plates can be modeled using a thermoviscoplastic constitutive model, e.g. a NortonHoff model, which characterizes the deviatoric component of the mechanical constitutive equation. Thus, only the thermal equation has to be solved for the tool. Using a NortonHoff constitutive model for the plates, the deviatoric component of the Cauchy stress tensor takes the form:

(84) 
where is the deviatoric component of the deformation rate, is the temperature dependent viscosity parameter, and is a temperature dependent rate sensitivity parameter. The mechanical dissipation rate per unit of spatial volume is computed as

(85) 
where is a parameter that represents the fraction of the total plastic dissipation rate converted into heat.
In the Solid Mechanics approach, two different constitutive models can be considered ([22], [23], [24]). A first option would be to consider a thermorigidplastic model, e.g. the NortonHoff model introduced above, to characterize the deviatoric component of the Cauchy stress tensor, and a thermoelastic model, written in rate form, to characterize the constitutive equation for the pressure and, thus, the volumetric component of the Cauchy stress tensor. Then, the pressure rate, under isothermal conditions, can be written as,

(86) 
where is the temperature dependent bulk modulus. Note that with this constitutive model, residual stresses can not be computed using a straightforward procedure. Alternatively, a thermoelastoviscoplastic constitutive model suitable for finite deformations can be used [94]. In this case, residual stresses can be computed following a straightforward procedure. On the other hand, plastic internal variables have to be integrated and the data transfer of the internal variables has to be performed.
In the Fluid Mechanics approach, the material is considered to be fully incompressible. The incompressibility constraint is incorporated into the equations to be solved, where the pressure plays the role of a Lagrange multiplier. Integration of the plastic strain tensor and data transfer of internal variables can be avoided. Residual stresses cannot be computed using a straightforward procedure. A twostep strategy to compute residual stresses within a Fluid Mechanics approach, based on a combination of local and global level analysis, is shown in Section 8.
Two different thermomechanical contact approaches are considered.
In the Solid Mechanics approach, fullstick thermomechanical contact conditions are considered between the tool and the workpieces. Thus, the displacement and temperature fields are continuous through the contact interface between the tool and the workpieces ([22], [23], [24]).
In the Fluid Mechanics approach, a thermofrictional contact model is considered ([22], [23], [24], [50]). Using a Norton thermofrictional contact model, the tangential component of the traction vector at the contact interface between the tool and the workpieces is given by,

(87) 
where and are the tangential component of the traction vector and the variation of the relative tangential velocity, respectively, at the contact interface between the tool and the workpieces, is the temperature dependent material consistency parameter, and is the strain rate sensitivity parameter. The heat flux generated by friction at the contact interface between the tool and the workpieces is split into a part absorbed by the tool, denoted as , and a part absorbed by the workpieces, denoted as , given, respectively, by ([22], [23], [24], [50])

(88) 
where the amount of heat absorbed by the tool, , and by the workpieces, , depends on the thermal diffusivity of the two materials in contact, and take the form

(89) 
where the thermal diffusivity parameter is given by

(90) 
Alternatively, as a limit case, fullstick thermomechanical contact conditions between the tool and the workpieces can be also considered. In this case the velocity and temperature fields are continuous through the contact interface between the tool and the workpieces.
In this section, two different formulations for the numerical simulation of FSW processes at local level have been compared. The first formulation, based on a Solid Mechanics approach, considers a quasiincompressible thermoelastoplastic material model (a thermorigidplastic deviatoric model as a particular case) and is written in terms of the displacement and temperature fields. The second one, based on a Fluid Mechanics approach, considers a fully incompressible thermorigidplastic material model and is written in terms of the velocity, pressure and temperature fields. The two approaches use different AKF and solution strategies. The solution scheme of the governing equations using an ALE formalism is also different for the two approaches. While remeshing and data transfer operations are avoided in the Fluid Mechanics approach, the Solid Mechanics approach requires to perform data transfer operations linked to the ALE formalism, as well as periodically remeshing and data transfer operations in the transition zone, to avoid excessive mesh distortion. Then, from a computational point of view, the Fluid Mechanics approach is more efficient than the Solid Mechanics one. One of the intrinsic drawbacks linked to the Fluid Mechanics approach is that residual stresses cannot computed using a straightforward procedure, while this is not the case if a Solid Mechanics approach, with a thermoelastoplastic consitutive model, is used.
Material flow visualization in horizontal and vertical directions. 
Figure 29: Material flow visualization in horizontal and vertical directions. 
Material points distribution. 
Figure 30: Material points distribution. 
In this method, firstly, a set of points representing the material points (tracers) are distributed in the domain and then, an Ordinary Differential Equation (ODE) for the computation of material displacement at a postprocess level must be solved (Figure particle distribution). Each particle position is integrated from:

(91) 
with the initial condition

Tracer position in time 
Figure 31: Tracer position in time 
where is the position of the material points at time and is the velocity of the tracer in the position and time . The tracer velocity is obtained interpolating the nodal velocity of the background mesh. To this end, firstly, a search algorithm is necessary in order to identify the element containing the tracer (Figure 31). To solve Eq. (91), there exist a large amount of integration techniques ranging from the simple first order Backward Euler (BE) scheme to higher order RungeKutta schemes. Among them three wellknown methods are chosen and compared: Backward Euler with Substepping (BES), Fourth order RungeKutta (RK4) and Back and Forth Error Compensation and Correction methods (BFECC). These methodologies, widely used in fluid dynamics, are suitable and robust tools to study the FSW problem allowing for a clear visualization of the material movement at the stirzone leading to a better understanding of the welding process.
According to the firstorder accurate BES method, the particle position at timestep is computed in substeps as:

(92) 
where is the timestep and is the substep. The velocity of the tracer at the current substep , , is computed as:

(93) 
Even if the accuracy of the method is enhanced by a substepping technique, it is still firstorder accurate
The steps of the BFECC method for a rotational velocity field 
Figure 32: The steps of the BFECC method for a rotational velocity field 
According to the thirdorder accurate, , BFECC method, the error arising from backward and forward advection of the particle position is used to modify the initial position. Denoting with the first order upwinding integration operator, BFECC consists of a prediction/correction procedure: the initial position of the tracer, , is advected forward to , and then backward to , being the backward advection operator. The difference between and the initial position defines the error, , of the time integration scheme. This error is used later to correct the starting position, , for the next advection step, achieving the final solution, . The steps are the following:

(94) 
where involves the computation of the velocity field of the tracer at its current position, interpolating from the nodal values (Figure 32).
Fourthorder RungeKutta method. In each step the derivative is evaluated four times: once at the initial point, twice at trial midpoints, and once at a trial endpoint. From these derivatives the final value (shown as a filled dot) is calculated 
Figure 33: Fourthorder RungeKutta method. In each step the derivative is evaluated four times: once at the initial point, twice at trial midpoints, and once at a trial endpoint. From these derivatives the final value (shown as a filled dot) is calculated 
Zalesask's benchmark after 10 revolutions. 
Figure 34: Zalesask's benchmark after 10 revolutions. 
Velocity effect on the weld quality. 
Figure 35: Velocity effect on the weld quality. 
Material particles after the weld with a threaded pin. 
Figure 36: Material particles after the weld with a threaded pin. 
According to the forthorder accurate RK4 method, the particle position at timestep is computed from the advection of the initial position by four weighted incremental displacement at intermediate timesteps (Figure 33):

(95) 
where the incremental displacements are computed as

(96) 
The RK4 method is the best choice among the other integration techniques because it can integrate exactly a circular trajectory, which is the most common particle path in FSW. Figure 34 shows the comparison between three represented methods in application to the Zalesask's problem. One can see that the RK4 method presents the exact result after 10 revolutions of the cutout disk. Figure 35 demonstrates the velocity effect on the final weld quality. It demonstrates that the rotational and advancing velocities can not be selected arbitrarily. Figure 36 shows the final state of the material particles (a colorful set) after FSW process.
Coupling strategy at local and global level. 
Figure 37: Coupling strategy at local and global level. 
The first step consists of computing the heat generated by the FSW tool at the thermomechanically affected zone (TMAZ). This is performed within the local level analysis.
In the local level analysis, the focus of the simulation is a small domain (typically, its size is the double of the pinshoulder) including the TMAZ. Because of this, the size of the local domain is not necessarily the whole component. It can be much smaller as big as the HAZ surrounded by artificial boundary (Figure 25). The local model provides a highly detailed description of the heat generation in order to capture accurately the temperature field. At local level, the analysis data are the pin geometry, rotation and advancing speeds, tool pressure and thermophysical material properties. The process phenomena studied are the relationship between the welding parameters, the contact mechanisms in terms of applied normal pressure and friction coefficient, the pin geometry, the material flow within the HAZ (or stirring zone), its size and, eventually, the corresponding consequences on the microstructure evolution. The local level analysis requires a very fine mesh in the TMAZ as in this region (basically a narrow area close to the pin) the strain rate gradient is extremely steep. This strain rate localization is mainly due to the highly nonlinear behavior of the constitutive models that characterize the material. The output of the local analysis, which is the heat power () to be used at global level, is calculated once the steadystate is achieved. It is computed taking into account the heat generated from both friction () and plastic dissipation () integrated over the TMAZ domain:

(97) 
where is the power density introduced into the system. Thus, performing the local level analysis provides the input data for the global study.
Schematic of a global model with a moving heat source 
Figure 38: Schematic of a global model with a moving heat source 
TMAZ schematic 
Figure 39: TMAZ schematic 
Schematic diagram of heat source model 
Figure 40: Schematic diagram of heat source model 
The numerical simulation at global level needs an adhoc procedure in order to apply the power energy to the elements representing the ThermoMechanically Affected Zone (TMAZ) (Figure 39). Therefore, a search algorithm is required to identify those elements at each time step of the simulation. The total power input () delivered within the timestep, , is distributed among the elements of the TMAZ volume during the current timestep. The upper surface of the volume is a rectangle with dimension of advancing speed () time step () pin diameter (). The diameter of the profile of this volume () decreases parabolically with depth in the workpiece thickness () (Figure 40):
