Local and global approaches to Friction Stir Welding processes

==Local and global approaches to Friction Stir Welding processes ==

'''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'''

==Abstract==

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 sub-domains is introduced for  the numerical modeling. As a result, the analysis can deal with complex  (non-cylindrical) pin-shapes and the extremely large deformation of the  material at the HAZ without requiring any remeshing or remapping tools.

A fully coupled thermo-mechanical 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 visco-plastic 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 thermo-elasto-  viscoplastic constitutive model. The analysis at local level is  characterized by a rigid thermo-visco-plastic constitutive model. Different  thermally coupled (non-Newtonian) fluid-like models as Norton-Hoff, Carreau  or Sheppard-Wright, among others are tested.

To better understand the material flow pattern in the stirring zone, a  (Lagrangian based) particle tracing is carried out while post-processing FSW  results.

A coupling strategy between the analysis of the process zone nearby the  pin-tool (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.

==1 Introduction==

===1.1 Industrial background of FSW===

Friction stir welding (FSW) is a solid state joining technology in which no  gross melting of the welded material takes place. It is a relatively new  technique (developed by The Welding Institute (TWI), in Cambridge, UK, in  1991) widely used over the past decades for joining aluminium alloys.  Recently, FSW has been applied to the joining of a wide variety of other  metals and alloys such as magnesium, titanium, steel and others. FSW is  considered to be the most significant development in metal joining in  decades and, in addition, is a "green" technology due to its energy  efficiency, environmental friendliness, and versatility. This process offers  a number of advantages over conventional joining processes (such as e.g.  fusion welding). The main advantages, often mentioned, include: (a) absence  of the need for expensive consumables such as a cover gas or flux; (b) ease  of automation of the machinery involved; (c) low distortion of the  work-piece; and (d) good mechanical properties of the resultant joint <span id='citeF-91'></span>[[#cite-91|[91]]], <span id='citeF-120'></span>[[#cite-120|[120]]]. Additionally, since it allows avoiding all  the problems associated to the cooling from the liquid phase, issues such as  porosity, solute redistribution, solidification cracking and liquation  cracking are not encountered during FSW<span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]]. In general, FSW has been found to produce a low  concentration of defects and is very tolerant to variations in parameters  and materials. Furthermore, since welding occurs by the deformation of  material at temperatures below the melting temperature, many problems  commonly associated with joining of dissimilar alloys can be avoided, thus  high-quality welds are produced. Due to this fact, it has been widely used  in different industrial applications where metallurgical characteristics  should be retained, such as in aeronautic, naval and automotive industry.   <div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-fsw-schematic-.png|600px|FSW process.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' FSW process.
|}

During FSW, the work-piece is placed on a backup plate and it is clamped  rigidly to eliminate any degrees of freedom (Figure [[#img-1|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 work-pieces, 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 <math display="inline">1-3   {{}^o}  </math> away from the direction of travel, although some tool designs allow it to  be positioned orthogonally to the work-piece (Figure [[#img-2|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 work-piece surface (see Figure [[#img-3|3]]). The  tool serves three primary functions, that is, heating of the work-piece,  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 work-piece top surface.<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-fig.png|600px|Tilt angle in FSW]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Tilt angle in FSW
|}
<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-4.png|600px|Schematic  representation of the friction stir welding process]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Schematic  representation of the friction stir welding process
|}

Depending on the geometrical configuration of the tool, material movement  around the pin can become complex, with severe gradients in temperature and  strain rate. This environment creates a challenge for modelers due to the  resulting coupled thermo-mechanical nature, the large deformation and strain  rates near the pin.<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-sch2.png|600px|Definition  of Friction Stir Welding zones Dong00]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-4|4]] that the work-pieces 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 <span id='citeF-27'></span>[[#cite-27|[27]]]. The periodic "onion flow" pattern  that is left behind as the tool advances is schematically illustrated in  Figure [[#img-4|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 post-weld properties depend on:

* Tool travel speed: influences total heat, porosity and weld quality.

* Tool rotation rate: influences total heat and weld quality.

* Tool design: shoulder diameter, scroll or concave shoulder, features  on the pin and pin length influence the extent of the material (Figure [[#img-5|5]]).

* Tool tilt: It influences the contact pressure. There exists lower  contact pressure (or incomplete contact) on the leading edge of the shoulder  due to tool tilt (typically between 0 <math display="inline">{{}^\circ }</math>   and 3 <math display="inline">{{}^\circ }</math>  ).

* Material thickness: influences cooling rate and through-thickness  temperature gradients.

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-pinshape.png|600px|Different pin shapes.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 key-point in FSW and its highly  non-linear 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 work-piece 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 work-piece around 80% to 90% of the melting temperature of the work-piece material, so that  welding defects are avoided <span id='citeF-33'></span>[[#cite-33|[33]]].

Moreover, analytical and numerical methods have a role to play, although  numerical methods dominate due to the accuracy and ease-to-use 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 two-dimensional (2-D) approximations to inherently three  dimensional (3-D) 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&#8211;-ideals that are rarely carried through in  practice.

<span id="fn-1"></span>
<span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]])  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.</span>

===1.2 Challenges for the simulation of FSW process===

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 set-ups 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. Computer-based 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  well-understood.

====1.2.1 Physical model====

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 work-pieces 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:

==Complex thermal behavior ==

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 work-pieces to the surrounding  environment while conduction heat flux occurs between the work-pieces and  the support.

==Non-linear behavior and localized nature ==

In FSW, the mechanical behavior is non-linear due to the high strain rates  and visco-plastic material. The strong non-linear region is limited to a  small area and the remaining part of the model is mostly linear. However,  the exact boundaries of the non-linear zone are not known a-priori. 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.

==Coupled nature ==

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,  temperature-dependent mechanical properties, temperature gradients in  work-pieces, etc.

An adequate physical model of the welding process must account for all these  phenomena including thermal, mechanical and coupling aspects.

====1.2.2 Numerical model====

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  thermo-mechanical 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 non-linear.  This has important implications upon the complexity of the numerical model.  Consequently, a robust and efficient numerical strategy is crucial for  solving such highly non-linear 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 visco-plastic 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 time-step 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.

==Mechanical problem==

A quasi-static 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 velocity-pressure  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 (Ladyzhenskaya-Babuska&#8211;Brezzi (LBB)  stability condition). Due to this, pressure instabilities appear if equal  velocity-pressure 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 thermo-mechanically coupled strain rate and temperature  dependent large deformation. According to the split of the stress tensor,  different rate-dependent constitutive models can be used for modeling of the  welding process. At typical welding temperatures, the large strain  deformation is mainly visco-plastic. Depending on the scope of the analysis,  rigid-visco-plastic or elasto-visco-plastic 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 non-linear  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.

==Thermal problem==

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 non-uniform 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 thermo-mechanical  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  thermo-mechanical 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.

==Kinematic framework==

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 history-dependent constitutive relations  are employed to describe elasto-plastic and visco-plastic 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  re-meshing. 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., fluid-fluid and fluid-solid  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.

==Coupled problem==

The numerical solution of the coupled thermo-mechanical problem involves the  transformation of an infinite dimensional transient system into a sequence  of discrete non-linear algebraic problems. This is achieved by means of the  FE spatial descritization procedure, a time-marching 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 time-stepping schemes, two types of strategies can be applied  to the solution of the coupled thermo-mechanical problems:

The first possibility is a monolithic (simultaneous) time-stepping 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 time-step as well as time-stepping  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 (block-iterative or  fractional-step), 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 non-linear 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 time-stepping  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.

==Particle tracing==

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 work-piece  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 time-step 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.

==Residual stresses==

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 high-quality 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 elasto-viscoplastic  constitutive model. Therefore, in this work, a local-global coupling  strategy is proposed in order to obtain the residual stress field.

===1.3 State of the art===

====1.3.1 Analytical and numerical modeling of FSW====

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 thermo-mehanical  modeling, global analysis and local modeling of the heat input and use of  different kinematic frameworks.

==Analytical modeling==

Thermal models are used to describe the heat flow in the FSW process.  Midling and Grong <span id='citeF-83'></span>[[#cite-83|[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 6082-T6 and AI-SiC, and indicated that a narrow  HAZ was obtained when a high specific power was used in conjunction with a  short duration heating cycle.

Grong <span id='citeF-66'></span>[[#cite-66|[66]]] described how Rosenthal's <span id='citeF-98'></span>[[#cite-98|[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 <span id='citeF-99'></span>[[#cite-99|[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. <span id='citeF-81'></span>[[#cite-81|[81]]], and Gould and Feng <span id='citeF-65'></span>[[#cite-65|[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. <span id='citeF-81'></span>[[#cite-81|[81]]] and was imposed on  a ring around the shoulder in Gould and Feng <span id='citeF-65'></span>[[#cite-65|[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 <span id='citeF-58'></span>[[#cite-58|[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 <math display="inline">{{}^\circ }</math>  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. <span id='citeF-112'></span>[[#cite-112|[112]]] also used an analytic thermal model, however the main emphasis  of this work was material flow.

Song and Kovacevic <span id='citeF-111'></span>[[#cite-111|[111]]] developed a mathematical model to describe  the detailed three-dimensional 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 non-uniform 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 <span id='citeF-57'></span>[[#cite-57|[57]]], a simple analytical  solution using series expansions for Cauchy momentum and energy equation-set  is obtained. The Power Law model takes into account the shear thinning fluid  and the Arrhenius-type relationship the temperature dependent viscosity.  Friction dissipation as an external heat source is considered.

==Numerical modeling==

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 work-piece, 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.

==Thermal modeling==

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 <span id='citeF-65'></span>[[#cite-65|[65]]] proposed a  simple heat transfer model to predict the temperature distribution in the  work-piece. Chao and Qi <span id='citeF-31'></span>[[#cite-31|[31]]], <span id='citeF-32'></span>[[#cite-32|[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 <span id='citeF-88'></span>[[#cite-88|[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 <span id='citeF-85'></span>[[#cite-85|[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. <span id='citeF-44'></span>[[#cite-44|[44]]], <span id='citeF-45'></span>[[#cite-45|[45]]] and  Frigaard et al. <span id='citeF-62'></span>[[#cite-62|[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/work-piece 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 <span id='citeF-84'></span>[[#cite-84|[84]]] and Russell and Sheercliff <span id='citeF-100'></span>[[#cite-100|[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. <span id='citeF-8'></span>[[#cite-8|[8]]]) were uncoupled or sequentially  coupled to the heat solvers, and limited by the computational power and  software capabilities of that time.

==Thermo-mechanical modeling==

More recently, a coupled thermomechanical modeling and simulation of the FSW  process can be found in Zhu and Chao <span id='citeF-122'></span>[[#cite-122|[122]]], Jorge Jr. and  Balancín <span id='citeF-74'></span>[[#cite-74|[74]]] or De Vuyst et al. <span id='citeF-49'></span>[[#cite-49|[49]]],   <span id='citeF-48'></span>[[#cite-48|[48]]]. Zhu et. al <span id='citeF-122'></span>[[#cite-122|[122]]] used a 3D  nonlinear thermal and thermo-mechanical 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  three-dimensional elasto-plastic thermo-mechanical 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. <span id='citeF-49'></span>[[#cite-49|[49]]],   <span id='citeF-48'></span>[[#cite-48|[48]]] used the coupled thermo-mechanical 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  Norton-Hoff rigid-visco-plastic 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 <span id='citeF-59'></span>[[#cite-59|[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 thermo-mechanical model into account. This  method permitted both transient and steady state analysis. The formulation  was developed in <span id='citeF-60'></span>[[#cite-60|[60]]] and <span id='citeF-67'></span>[[#cite-67|[67]]] to simulate the  different stages of the FSW process. Assidi <span id='citeF-10'></span>[[#cite-10|[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. <span id='citeF-55'></span>[[#cite-55|[55]]] and Chao et al. <span id='citeF-33'></span>[[#cite-33|[33]]].

In <span id='citeF-33'></span>[[#cite-33|[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 work-piece. To validate the  result, the temperature evolution was recorded in the tool and in the  work-piece. 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 work-piece the code WELDSIM was used.  It was a transient, nonlinear, 3D FE code. In this model only half of the  work-piece 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  work-piece 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 <span id='citeF-93'></span>[[#cite-93|[93]]], the Levenberg-Marquardt (LM) method is used  in order to perform a non-linear 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.

==Global level modeling==

Most of the above-mentioned 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 <span id='citeF-108'></span>[[#cite-108|[108]]], <span id='citeF-109'></span>[[#cite-109|[109]]], <span id='citeF-110'></span>[[#cite-110|[110]]], <span id='citeF-31'></span>[[#cite-31|[31]]], [Chao-Qi
1999], <span id='citeF-44'></span>[[#cite-44|[44]]]. Another technique involved estimating  the power input analytically and then artificially limiting the peak  temperature <span id='citeF-108'></span>[[#cite-108|[108]]], <span id='citeF-62'></span>[[#cite-62|[62]]] or introducing latent  heat effects <span id='citeF-108'></span>[[#cite-108|[108]]], <span id='citeF-109'></span>[[#cite-109|[109]]] to avoid over-predicting the  weld temperature. The most satisfactory approach involved measuring the weld  power experimentally and using this as an input to the model <span id='citeF-79'></span>[[#cite-79|[79]]], <span id='citeF-107'></span>[[#cite-107|[107]]], <span id='citeF-75'></span>[[#cite-75|[75]]], <span id='citeF-76'></span>[[#cite-76|[76]]], <span id='citeF-77'></span>[[#cite-77|[77]]].

Khandkar et al. <span id='citeF-76'></span>[[#cite-76|[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  pin-tool was correlated to input torque data obtained from experimental  investigation of butt-welding. The moving heat source included heat  generation due to torques at the interface between the tool shoulder and the  work-piece, the horizontal interface between the pin bottom and the  work-piece, and the vertical interface between the cylindrical pin surface  and the work-piece. Temperature-dependent properties of the weld-material  were used for the numerical modeling.

In Khandkar et al. <span id='citeF-77'></span>[[#cite-77|[77]]] and Hamilton et al. <span id='citeF-68'></span>[[#cite-68|[68]]] a  torque-based 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 energy-slip 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 <span id='citeF-78'></span>[[#cite-78|[78]]].Khandkar  et al. <span id='citeF-78'></span>[[#cite-78|[78]]] used a coupled thermo-mechanical 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 <span id='citeF-97'></span>[[#cite-97|[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, thermo-physical 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,  two-dimensional fluid dynamics based model that was used to make parametric  studies of variations in properties of the material to be welded (mechanical  and thermo-physical) and variations in welding parameters. This was done by  a non slip boundary condition at the tool work-piece interface. The  deformation behavior was based on deviatoric flow stress using the  Zener-Hollomon 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 thermo-mechanical models for the prediction of  the distribution of the residual stresses in the process of friction stir  welding. A steady-state simulation of FSW was carried out by Bastier et al.   <span id='citeF-12'></span>[[#cite-12|[12]]]. The simulation included two main steps. The first one  uses an Eulerian description of the thermo-mechanical problem together with  a steady-state algorithm detailed in <span id='citeF-47'></span>[[#cite-47|[47]]], in order to avoid  remeshing due to the pin motion. In the second step, a steady-state  algorithm based on an elasto-visco-plastic 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. <span id='citeF-82'></span>[[#cite-82|[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. <span id='citeF-92'></span>[[#cite-92|[92]]] used a numerical-experimental 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  out-of-plane displacements as input nodal boundary conditions.

Yan et al. <span id='citeF-118'></span>[[#cite-118|[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 <span id='citeF-61'></span>[[#cite-61|[61]]] used the cut-compliance and the inverse  weight-function methodologies for skin stringer FSW geometries via finite  element analysis to measure residual stresses.

Rahmati Darvazi and Iranmanesh <span id='citeF-96'></span>[[#cite-96|[96]]] presented a  thermo-mechanical model to predict the longitudinal residual stress applying  a so-called advancing retreating factor. The uncoupled thermo-mechanically  equations were solved using ABAQUS.

==Local level modeling==

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 visco-plastic 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 <span id='citeF-80'></span>[[#cite-80|[80]]], <span id='citeF-64'></span>[[#cite-64|[64]]] and <span id='citeF-63'></span>[[#cite-63|[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 <span id='citeF-80'></span>[[#cite-80|[80]]] (FEM with a Lagrangian mesh),  the material was assumed to be visco-plastic, with temperature-dependent  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 <span id='citeF-15'></span>[[#cite-15|[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 Eulerian-Lagrangian mesh.

A thermo-mechanical model for FSW was proposed by Dong et al. <span id='citeF-54'></span>[[#cite-54|[54]]].  There, an axis-symmetrical FE Lagrangian formulation was used. Ulysse in   <span id='citeF-113'></span>[[#cite-113|[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 stir-welding  process using 3D visco-plastic 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. <span id='citeF-9'></span>[[#cite-9|[9]]] used the CTH finite volume hydrocode coupled  to an advection-diffusion 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. <span id='citeF-115'></span>[[#cite-115|[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 <span id='citeF-116'></span>[[#cite-116|[116]]], <span id='citeF-117'></span>[[#cite-117|[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 elasto-plastic deformations and temperature  dependent material properties. However, the authors did not perform a fully  coupled thermo-mechanical 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 <span id='citeF-104'></span>[[#cite-104|[104]]] also used the CFD  commercial software FLUENT to model the 2D steady-state flow around a  cylindrical tool. The paper describes the progressive development of a  finite volume model in FLUENT that used a visco-plastic material whose  viscous properties were based on the Sellars-Tegart relationship. Even  though the model was 2D, heat generation and conduction were included. To  avoid over-predicting 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. <span id='citeF-13'></span>[[#cite-13|[13]]], <span id='citeF-14'></span>[[#cite-14|[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 Eulerian-Lagrangian mesh.

Schmidt and Hattel <span id='citeF-103'></span>[[#cite-103|[103]]] presented a 3D fully coupled  thermo-mechanical FE model in ABAQUS/Explicit using the ALE formulation and  the Johnson-Cook material law. The flexibility of the FSW machine was taken  into account connecting the rigid tool to a spring. The work-piece was  modeled as a cylindrical volume with inlet and outlet boundary conditions. A  rigid back-plate 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 steady-state (14 days on a 3 GHz Pentium PC to reach only 10  seconds of model time).

The model developed by Chen and Kovacevic in <span id='citeF-34'></span>[[#cite-34|[34]]] uses the  commercial FEM software ANSYS for a thermo-mechanically coupled Lagrangian  finite element model of aluminium alloy AA-6061-T6. 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 work-piece material.  Finally the model predicted the welding forces in the x, y and z directions.

Nikiforakis <span id='citeF-89'></span>[[#cite-89|[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 thermo-fluid analysis was performed. The rotation of the  tool was handled through the use of the overlapping grid method. A  rigid-visco-plastic material law was used and sticking contact at the tool  work-piece interface was assumed. Hence, heating was due to plastic  deformation only.

Heurtier et al. <span id='citeF-69'></span>[[#cite-69|[69]]] used a 3D semi-analytical coupled  thermo-mechanical 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 micro-hardness in various weld zones are computed and  compared to experimental results obtained on a AA 2024-T351 alloy FSW joint.

==Kinematic framework==

Among the distinctive local level studies, Cho et al. <span id='citeF-43'></span>[[#cite-43|[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 re-meshing was employed in <span id='citeF-35'></span>[[#cite-35|[35]]] while similar approaches were applied in <span id='citeF-19'></span>[[#cite-19|[19]]] and <span id='citeF-20'></span>[[#cite-20|[20]]], which are not numerically efficient.

Buffa et al. <span id='citeF-19'></span>[[#cite-19|[19]]] using the commercial FE software DEFORM-3D,  proposed a 3D Lagrangian, implicit, coupled thermo-mechanical numerical  model for the simulation of FSW processes, using a rigid-visco-plastic  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 <span id='citeF-86'></span>[[#cite-86|[86]]] and <span id='citeF-87'></span>[[#cite-87|[87]]] employed a control  volume approach for discretization of the FSW domain. They investigated  visco-plastic 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 non-Newtonian 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 (non-fluid) solvers, which only predict weld temperature.

Recently, Assidi et al. in <span id='citeF-11'></span>[[#cite-11|[11]]] used an ALE formulation  implemented into the Forge3 software with a splitting approach  and an adaptive re-meshing scheme based on error estimation. In <span id='citeF-21'></span>[[#cite-21|[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. <span id='citeF-21'></span>[[#cite-21|[21]]] simulated the welding process using a continuous  rigid-viscoplastic finite element model in a single block approach through  the Lagrangian implicit software, DEFORM-3DTM. Then, the temperature  histories extracted at each node of the model were transferred to another  finite element model considering elasto-plastic 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. <span id='citeF-101'></span>[[#cite-101|[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/work-piece interface.

The more recent works performed in this direction are those presented in  this paper. Agelet de Saracibar et al. (<span id='citeF-3'></span>[[#cite-3|[3]]], <span id='citeF-4'></span>[[#cite-4|[4]]]), Agelet de Saracibar et al. (<span id='citeF-5'></span>[[#cite-5|[5]]], <span id='citeF-6'></span>[[#cite-6|[6]]]), Chiumenti et al. <span id='citeF-40'></span>[[#cite-40|[40]]], and Dialami et al. ( <span id='citeF-50'></span>[[#cite-50|[50]]], <span id='citeF-51'></span>[[#cite-51|[51]]], <span id='citeF-52'></span>[[#cite-52|[52]]]) used a sub-grid  scale finite element stabilized mixed velocity/pressure/temperature  formulation for coupled thermo-rigid-plastic models, using Eulerian and  Arbitrary Lagrangian Eulerian (ALE) formalisms, for the numerical simulation  of FSW processes. They used ASGS and OSGS methods and quasi-static sub-grid  scales, neglecting the sub-grid scale pressure and using the finite element  component of the velocity in the convective term of the energy balance  equation. Chiumenti et al. <span id='citeF-40'></span>[[#cite-40|[40]]], Dialami et al. <span id='citeF-50'></span>[[#cite-50|[50]]], and Chiumenti et al. <span id='citeF-41'></span>[[#cite-41|[41]]] developed an  apropos kinematic framework for the numerical simulation of FSW processes  and compared with a solid approach in <span id='citeF-22'></span>[[#cite-22|[22]]] and <span id='citeF-24'></span>[[#cite-24|[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 <span id='citeF-30'></span>[[#cite-30|[30]]] developed by the  authors at the International Centre for Numerical Methods in Engineering  (CIMNE). Chiumenti et al. <span id='citeF-41'></span>[[#cite-41|[41]]] used a novel stress-accurate  FE technology for highly nonlinear analysis with incompressibility  constraints typically found in the numerical simulation of FSW processes.  They used a mixed linear piece-wise 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.

===1.4 Outline===

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 thermo-elasto-visco-plastic (global level  analysis) or thermo-rigid-visco-plastic (local level analysis).

Section 4 describes the discrete FE modeling and sub-grid scale  stabilization to deal with incompressibility and heat convection dominated  problems. The multiscale stabilization method is introduced and an  approximation of the sub-grid scale variables together with the  stabilization parameters is given. Algebraic Sub-grid Scale (ASGS) and  Orthogonal Sub-grid 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 Arbitrary-Lagrangian-Eulerian kinematic framework, the  overall computational domain is divided into sub-domains 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 post-process level.

Section 8 is devoted to the description of the Local-global 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.

==2 Thermal problem==

The governing equation representing the thermal problem is the ''balance  of energy'' equation. This equation controls the temperature evolution and  can be stated as (<span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-28'></span>[[#cite-28|[28]]],   <span id='citeF-38'></span>[[#cite-38|[38]]]):

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _{o}c\frac{dT}{dt}=\dot{R}-\mathbf{\nabla }\cdot \mathbf{q}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\rho _{o},</math> <math display="inline">c,T</math>, <math display="inline">\dot{R}</math> and <math display="inline">\mathbf{q}</math> 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, <math display="inline">\dot{R}</math>:

* Power dissipated through plastic deformation (i.e. local analysis ).

* Known input (i.e. global analysis ).

===2.1 Thermal problem at the global level===

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-global-temp.png|600px|Temperature field in a FSW  process (Global level)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Temperature field in a FSW  process (Global level)
|}
<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-temp-glob-line.png|600px|Temperature evolution obtained from global level  analysis compared with experiment]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Temperature evolution obtained from global level  analysis compared with experiment
|}
<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-fsw-local.png|600px|Temperature contour field in FSW process (Local level).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' Temperature contour field in FSW process (Local level).
|}
<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-temp-loc-line.png|600px|Temperature evolution obtained from local level analysis  compared with experiment]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _{o}c\dfrac{\partial T}{\partial t}=\dot{R}-\mathbf{\nabla }\cdot   \mathbf{q}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

In global level analysis, the heat source is the (known) power input for the  problem, obtained by solving the local level analysis.

Figure [[#img-6|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 [[#img-7|7]].

===2.2 Thermal problem at the local level===

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

<span id="eq-3"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\rho _{o}c\left(\dfrac{\partial T}{\partial t}+\mathbf{v}\cdot \mathbf{ \nabla }T\right)=\dot{R}-\mathbf{\nabla }\cdot \mathbf{q}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

where <math display="inline">\mathbf{v}</math> 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, <math display="inline">\mathbf{\dot{e}}</math>, and the deviatoric stresses, <math display="inline">\mathbf{s}</math>,  as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{R}=\gamma \mathbf{s:\dot{e}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

where <math display="inline">\gamma \simeq 90%</math> 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 [[#img-8|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 [[#img-9|9]].

===2.3 Boundary conditions===

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-bc.png|600px|Thermal boundary condition]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' Thermal boundary condition
|}
Let us denote by <math display="inline">\Omega </math> an open  and bounded domain. The boundary <math display="inline">\partial \Omega </math> can be split into <math display="inline"> \partial \Omega _{q}</math> and <math display="inline">\partial \Omega _{T}</math> such that <math display="inline">\partial \Omega  =\partial \Omega _{q}\cup </math> <math display="inline">\partial \Omega _{T}</math>, where fluxes (on <math display="inline"> \partial \Omega _{q}</math>) and temperatures (on <math display="inline">\partial \Omega _{T}</math>) are  prescribed for the heat transfer analysis as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{ccc}\mathbf{q\cdot n}=\bar{q} & \hbox{ on} & \partial \Omega _{q} \\   T=\bar{T} & \hbox{ on} & \partial \Omega _{T} \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

In figure [[#img-10|10]], <math display="inline">\Omega _{1}</math> and <math display="inline">\Omega _{2}</math> represent the  tool and the work-piece domains, respectively.

The initial condition for the transient thermal problem in terms of the  initial temperature field is: <math display="inline">T\left(t=0\right)=T_{o}</math>.

On free surfaces, the heat flux is dissipated through the boundaries to the  environment by heat convection, expressed by Newton's law as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{conv}=h_{conv}(T-T_{env})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

where <math display="inline">h_{conv}</math> is the heat transfer coefficient by convection, <math display="inline">T_{env}</math>  is the surrounding environment temperature and <math display="inline">T</math> is the temperature of the  body surface.

Another heat dissipation mechanism is heat loss due to radiation. Heat  radiation flux is computed using the Stefan-Boltzmann law:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{rad}=\sigma _{0}\varepsilon (T^{4}-T_{env}^{4})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

where <math display="inline">\sigma _{0}</math> is the Stefan&#8211;Boltzmann constant and <math display="inline">\varepsilon </math> is  the emissivity factor.

===2.4 Thermal constitutive model===

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, <math display="inline">\nabla T</math>, and the  (temperature dependent) thermal conductivity, <math display="inline">k</math>, as:

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{q}=-k\,\nabla T    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

===2.5 Friction model===

The thermal exchanges at the contact boundary (<math display="inline">\partial \Omega _{c}</math> in  figure [[#img-10|10]]) can also result from a friction type dissipation  process (Figure [[#img-11|11]]). The heat generated by frictional  dissipation at the contact interface is absorbed by the two bodies in  contact according to their respective thermal diffusivity. In this section  the frictional contact is described by two models: Coulomb's law used when  Lagrangian setting is considered in global studies; and Norton's law  considered when Eulerian/ALE setting is used in local studies.<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-stick2.png|600px|Velocity (top) and temperature fields (bottom)  obtained with fully stick condition between pin and work-piece.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 11:''' Velocity (top) and temperature fields (bottom)  obtained with fully stick condition between pin and work-piece.
|}

====2.5.1 Coulomb's friction law====

Adopting the classical friction model based on ''Coulomb's law'', the so  called ''slip function'' is defined as <span id='citeF-7'></span>[[#cite-7|[7]]]:

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\phi \left(\mathbf{t}_{N},\mathbf{t}_{T}\right)=\left\Vert \mathbf{t}_{T}\right\Vert -\eta \left\Vert \mathbf{t}_{N}\right\Vert \leq 0    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

where <math display="inline">\eta \left(T,\Delta \mathbf{u}_{c}\right)</math> is the friction  coefficient, which can result in a nonlinear function of temperature and  slip displacement <math display="inline">\Delta \mathbf{u}_{c}</math>. <math display="inline">\mathbf{t}_{N}</math> and <math display="inline">\mathbf{t}_{T}</math> are the normal and tangential components of the traction vector, <math display="inline"> \mathbf{t}_{c}=\mathbf{\sigma }\cdot \mathbf{n}</math>, at the contact interface,  respectively:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{t}_{N} =\left(\mathbf{n}\otimes \mathbf{n}\right)\cdot \mathbf{t}_{c}=\left(\mathbf{t}_{c}\cdot \mathbf{n}\right)\mathbf{n} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|-
| style="text-align: center;" | <math> \mathbf{t}_{T} =\left(\mathbf{I}-\mathbf{n}\otimes \mathbf{n}\right) \cdot \mathbf{t}_{c}=\mathbf{t}_{c}-\mathbf{t}_{N}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}
|}

where <math display="inline">\mathbf{n}</math> is the unit vector normal to the contact interface.

This given, both stick and slip mechanisms can be recovered using the  unified format:

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\Vert \mathbf{t}_{T}\right\Vert =\varepsilon _{T}\left(\left\Vert  \Delta \mathbf{u}_{T}\right\Vert -\gamma \frac{\partial \phi }{\partial  \left\Vert \mathbf{t}_{T}\right\Vert }\right)=\varepsilon _{T}\left( \left\Vert \Delta \mathbf{u}_{T}\right\Vert -\gamma \right)   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

together with the Kuhn-Tucker conditions defined in terms of the slip  function, <math display="inline">\phi </math>, and the slip multiplier, <math display="inline">\gamma </math>, as:

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{c}\phi \leq 0 \\   \gamma \geq 0 \\   \phi \dot{\gamma }=0 \end{array}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where <math display="inline">\varepsilon _{T}</math> is a penalty parameter (regularization of the  Heaviside function).

<math display="inline">\Delta \mathbf{u}_{T}</math> and <math display="inline">\Delta \mathbf{u}_{N}</math> are the tangential and  the normal components of the total relative displacement, <math display="inline">\Delta \mathbf{u}_{c},</math> computed as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Delta \mathbf{u}_{T} =\left(\mathbf{I}-\mathbf{n}\otimes \mathbf{n}\right)\cdot \Delta \mathbf{u}_{c} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|-
| style="text-align: center;" | <math> \Delta \mathbf{u}_{N} =\Delta \mathbf{u}_{c}-\Delta \mathbf{u}_{T}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}
|}

* The ''sliding condition'' allows for relative slip between the  contact surfaces. The tangential traction vector, <math display="inline">\mathbf{t}_{T}</math>, is  computed from ([[#eq-9|9]]) assuming <math display="inline">\phi =0</math>, as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{t}_{T}=\eta \left\Vert \mathbf{t}_{N}\right\Vert \frac{\Delta   \mathbf{u}_{T}}{\left\Vert \Delta \mathbf{u}_{T}\right\Vert }  

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

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

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{t}_{N}=\varepsilon _{N}\Delta \mathbf{u}_{N}    

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

where <math display="inline">\varepsilon _{N}</math> is the normal penalty parameter (not necessarily  equal to <math display="inline">\varepsilon _{T}</math>), which enforces non-penetration in the normal  direction.

* The ''stick condition'' is obtained for <math display="inline">\dot{\gamma }=0</math>, when the  contact surfaces are sticked to each other and there is no relative slip  between them. The tangential traction vector reads:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{t}_{T}=\varepsilon _{T}\Delta \mathbf{u}_{T}  

</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

while the normal component of the traction vector is given by Eq. ([[#eq-17|17]]).

====2.5.2 Norton's friction law====

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, <math display="inline"> \mathbf{v}</math>, it is very convenient to use a Norton type friction law.

The tangential component of the traction vector at the contact interface, <math display="inline"> \mathbf{t}_{T}</math>, is defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{t}_{T}=\eta _{eq}\mathbf{n}_{T}=a\left(T\right)\left\Vert \Delta   \mathbf{v}_{T}\right\Vert ^{q}\mathbf{n}_{T}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

where <math display="inline">\eta _{eq}\left(T,\Delta \mathbf{v}_{T}\right)=a\left(T\right) \left\Vert \Delta \mathbf{v}_{T}\right\Vert ^{q}</math> is the equivalent  friction coefficient. <math display="inline">a\left(T\right)</math> is the (temperature dependent)  material consistency, <math display="inline">0\leq q\leq 1</math> is the strain rate sensitivity and <math display="inline"> \mathbf{n}_{T}=\dfrac{\Delta \mathbf{v}_{T}}{\left\Vert \Delta \mathbf{v}_{T}\right\Vert }</math> 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>q_{frict}^{\left(1\right)} =-\varsigma ^{\left(1\right)}\mathbf{t}_{T}\cdot \Delta \mathbf{v}_{T}=-\varsigma ^{\left(1\right)}a\left( T\right)\left\Vert \Delta \mathbf{v}_{T}\right\Vert ^{q+1} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|-
| style="text-align: center;" | <math> q_{frict}^{\left(2\right)} =-\varsigma ^{\left(2\right)}\mathbf{t}_{T}\cdot \Delta \mathbf{v}_{T}=-\varsigma ^{\left(2\right)}a\left( T\right)\left\Vert \Delta \mathbf{v}_{T}\right\Vert ^{q+1}   </math>
|}
|}

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, <math display="inline">\varsigma ^{\left(1\right)}</math>, and by the second body, <math display="inline"> \varsigma ^{\left(2\right)}</math>, depends on the thermal diffusivity, <math display="inline">\alpha = \dfrac{k}{\rho _{0}c}</math>, of the two materials in contact as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\varsigma ^{\left(1\right)} =\dfrac{\alpha ^{\left(1\right)}}{\alpha  ^{\left(1\right)}+\alpha ^{\left(2\right)}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|-
| style="text-align: center;" | <math> \varsigma ^{\left(2\right)} =\dfrac{\alpha ^{\left(2\right)}}{\alpha  ^{\left(2\right)}+\alpha ^{\left(1\right)}}   </math>
|}
|}

The more diffusive is the material of one part in comparison with the other  part, the more heat is absorbed by it.

===2.6 Weak form of the thermal problem===

Considering Eqs. ([[#eq-2|2]]) and ([[#eq-8|8]]), the weak form of  the thermal problem at global level is defined over the integration domain <math display="inline"> \Omega </math> and its boundary <math display="inline">\partial \Omega </math> as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int \limits _{\Omega }\left(\left[\rho _{o}c\dfrac{\partial T}{\partial t} \right] \delta T\right)dV+ </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|-
| style="text-align: center;" | <math> \int \limits _{\Omega }\left[k\nabla T\cdot \nabla \left(\delta T\right)  \right] dV =W_{ther}\hbox{     }\forall \delta T   </math>
|}
|}

where <math display="inline">\delta T</math> is the test function of the temperature field, while the  thermal work, <math display="inline">W_{ther}</math>, is defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>W_{ther}=\int \limits _{\Omega }\left(\dot{R}\mathbf{ }\delta T\right) dV-\int \limits _{\partial \Omega q}\left(\bar{q}\delta T\right) dS-\int \limits _{\partial \Omega _{c}}\left(q_{frict}\delta T\right)dS  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|}

Note that the term due to friction (<math display="inline">q_{frict}</math>) appears only at local level.

Similarly, considering Eqs. ([[#eq-3|3]]) and ([[#eq-8|8]]), the  weak form of the problem at local level is written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\int \limits _{\Omega }\left(\left[\rho _{o}c\left(\dfrac{\partial T}{ \partial t}+\mathbf{v}\cdot \mathbf{\nabla }T\right)\right] \delta  T\right)dV+ </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|-
| style="text-align: center;" | <math> \int \limits _{\Omega }\left[k\nabla T\cdot \nabla \left(\delta T\right)  \right] dV =W_{ther}\hbox{     }\forall \delta T   </math>
|}
|}

===2.7 Discrete weak form of the thermal problem===

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

* Global form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\int \limits _{\Omega }\left[\left(\rho _{o}c\dfrac{\partial T_{h}}{\partial  t}\right) \delta T_{h}\right] dV </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|-
| style="text-align: center;" | <math> +\int \limits _{\Omega } k\nabla T_{h}\mathbf{\cdot \nabla }\left(\delta  T_{h}\right)dV  =W_{ther}^{ext}\left(\delta T_{h}\right)\hbox{     }\forall \delta T_{h}   

</math>
|}
|}

* Local form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\int \limits _{\Omega }\left[\left(\rho _{o}c\left(\dfrac{\partial T_{h}}{ \partial t}+\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right)\right) \delta T_{h}\right] dV </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|-
| style="text-align: center;" | <math> +\int \limits _{\Omega } k\nabla T_{h}\mathbf{\cdot \nabla }\left(\delta  T_{h}\right)dV  =W_{ther}^{ext}\left(\delta T_{h}\right)\hbox{     }\forall \delta T_{h}   

</math>
|}
|}

where <math display="inline">T_{h}</math> is the discrete temperature field.

==3 Mechanical problem==

In this section, the general framework for the description of the mechanical  problem in FSW processes is presented. This framework is developed both for  local and global analysis. In global level analyses, the effect of a moving  heat source on the entire structure is studied. The heat source is  introduced into the system and moves along the weld line while the structure  is fixed. In this case the Lagrangian framework is used for the definition  of the problem, as the reference frame is attached to the structure. The  moving heat source generates thermal deformation in the structure.  Therefore, definition of the problem in terms of displacements is a natural  choice for the computation of the mechanical problem. On the other side, in  local level analysis, the focus of the study is the heat affected zone. The  heat source is fixed together with the reference frame while the structure  moves with the imposed velocity (relative movement). In this case, the  Eulerian model is a suitable choice and the problem can be conveniently  defined in terms of the velocity field.<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-mech-res.png|600px|Mechanical results in a FSWprocess (local analysis).]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 12:''' Mechanical results in a FSWprocess (local analysis).
|}

In both cases, assuming quasi-static conditions, the local form of the '' balance of momentum equation, ''also known as ''Cauchy's equation of  motion,'' is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\nabla }\cdot \mathbf{\sigma }+\rho _{o}\mathbf{b}=0  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

where <math display="inline">\mathbf{\sigma }</math> is the Cauchy's stress tensor, <math display="inline">\mathbf{b}</math> is the  body forces vector per unit of mass, <math display="inline">\rho _{o}</math> is the density in the  reference configuration and <math display="inline">\mathbf{\nabla }\cdot \left(\circ \right)</math> 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, <math display="inline">\mathbf{\sigma }</math>, into volumetric and  deviatoric parts:

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\sigma }=p\mathbf{I}+\mathbf{s}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

where <math display="inline">p=\frac{1}{3}\mathrm{tr}\left(\mathbf{\sigma }\right)</math> is the  pressure and <math display="inline">\mathbf{s}=\mathrm{dev}\left(\mathbf{\sigma }\right)</math> is the  stress deviator. Similarly, the strain tensor, <math display="inline">\mathbf{\varepsilon ,}</math> can  be split into volumetric, <math display="inline">e_{vol},</math> and the deviatoric, <math display="inline">\mathbf{e,}</math>  parts.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\varepsilon }=\frac{1}{3}e_{vol}\mathbf{I}+\mathbf{e}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

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

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\nabla }\cdot \mathbf{s+\nabla }p+\rho _{o}\mathbf{b}=\mathbf{0}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

where <math display="inline">p</math> and <math display="inline">\mathbf{s}</math> are defined according to the constitutive  equations selected for the local or global analysis, respectively.

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-sl.png|600px|Velocity streamlines around a  triangular pin (left) and a trivex pin (right)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 13:''' Velocity streamlines around a  triangular pin (left) and a trivex pin (right)
|}

<div id='img-14'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-plastic-dissipation.png|600px|Plastic dissipation under a  FSWtool.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 14:''' Plastic dissipation under a  FSWtool.
|}

<div id='img-15'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-residu_0.png|600px|Longitudinal stress (XX) field  in a global level simulation of FSW after cooling.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 15:''' Longitudinal stress (XX) field  in a global level simulation of FSW after cooling.
|}
Figure [[#img-12|12]] shows  the results of the mechanical simulation at the local level. Pressure,  velocity, strain-rate 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 [[#img-13|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 [[#img-14|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 [[#img-15|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.

===3.1 Constitutive model===

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.

====3.1.1 Mechanical constitutive model at global 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, <math display="inline">\mathbf{s}\left(\mathbf{u},T\right)</math> , is generally defined through the constitutive equation as a function of  displacements <math display="inline">\mathbf{u}</math> and temperature <math display="inline">T.</math> 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 elasto-plastic (at room temperature) to pure  viscous (close to the melting point). This evolution can be defined by a  thermo-elasto-visco-plastic constitutive model. Close to the melting point,  the viscous behavior dominates and the elastic limit gradually decreases.

To this end, it can be written:

<span id="eq-31"></span>
<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>p =K\,e^{e}=K\,\,\left(e\mathbf{-}e^{\theta \,\,}\right)=K\,\,\left( \nabla \cdot \mathbf{u}-e^{\theta \,\,}\right)  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|-
| style="text-align: center;" | <math> \mathbf{s} =G\mathbf{e}^{e}=\,G\,\left(\,\,\mathbf{\,e-e}^{vp}\right)   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}
|}

where <math display="inline">K\left(T\right)</math> is the bulk modulus, which controls the material  compressibility and <math display="inline">G</math> is the shear modulus. Deformations <math display="inline">e^{\theta  \,\,},e^{e},\mathbf{e}^{vp}</math> and <math display="inline">\mathbf{e}^{e}</math> are the volumetric thermal  deformation, the elastic volumetric strain, the viscoplastic and the elastic  deviatoric strain, respectively. The thermal deformation <math display="inline">e^{\theta \,\,}</math>  is expressed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>e^{\theta \,\,}=3\left[\alpha \left(T\right)\left(T-T_{ref}\right) -\alpha \left(T_{0}\right)\left(T_{0}-T_{ref}\right)\right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

where <math display="inline">\alpha \left(T\right)</math> is the thermal expansion coefficient.

<div id='img-16'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-property.png|600px|Temperature dependent material  properties of 304L stainless steel.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 16:''' Temperature dependent material  properties of 304L stainless steel.
|}

Figure [[#img-16|16]] shows the temperature-dependent material property for  304L stainless steel. The mechanical material parameters can be obtained by  fitting the constitutive model with uniaxial, tensile stress-strain curves  obtained at different temperatures.

The viscoplastic deformation <math display="inline">\mathbf{\dot{e}}^{vp}</math> is defined through  appropriate evolution laws

<span id="eq-34"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\dot{e}}^{vp} =\gamma ^{vp}\,\dfrac{\mathbf{\partial }\Phi }{ \partial \mathbf{s}}=\gamma ^{vp}\mathbf{n}   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|-
| style="text-align: center;" | <math> \dot{\xi } =\gamma ^{vp}\,\sqrt{\dfrac{2}{3}\,}  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}
|}

where <math display="inline">\Phi \left(\mathbf{s,\,}T\right)</math> is the viscoplastic yield  surface, <math display="inline">\mathbf{n}=\dfrac{\mathbf{s}}{\left\Vert \mathbf{s}\right\Vert }</math>  is the unit normal to the yield surface and <math display="inline">\gamma ^{vp}\,=\left\langle   \dfrac{\Phi }{\eta }\right\rangle </math> is the visco-plastic parameter, when a  rate-dependent evolution law is assumed. <math display="inline">\eta \left(T\right)</math> is the  plastic viscosity associated to the visco-plastic model.

As in this model the elastic effects are not neglected, the residual  stresses can be evaluated.

Note that in the liquid-like phase, the visco-plastic potential

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Phi =\left\Vert \mathbf{s}\right\Vert -\tau \left\Vert \mathbf{e}^{vp}\right\Vert  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

degenerates into a purely-viscous potential <math display="inline">\Phi =\left\Vert \mathbf{s}\right\Vert </math>, when the coefficient <math display="inline">\tau \rightarrow 0</math>. As a consequence,  considering the definition of the visco-plastic multiplier, it transforms  into <math display="inline">\left\Vert \mathbf{s}\right\Vert =\gamma ^{vp}\eta </math>.  This means  that the visco-plastic model reduces to a Norton-Hoff model which is a  viscous relationship between the deviator of the stress tensor, <math display="inline">\mathbf{s}</math> ''',''' and the visco-plastic strain rate, <math display="inline">\mathbf{\dot{e}}^{vp}</math>,  governed by the viscous parameter, <math display="inline">\eta </math>.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{s}=\gamma ^{vp}\eta \hbox{ }\mathbf{n=}\eta \hbox{ }\mathbf{\dot{e}}^{vp}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

====3.1.2 Mechanical constitutive model at local level====

When the heat affected zone is the focus of the study, the local level  analysis in Eulerian format is preferred. The stress deviator <math display="inline">\mathbf{s}\left(\mathbf{v},T\right)</math> is defined as a function of temperature, <math display="inline">T,</math>  and the velocities, <math display="inline">\mathbf{v,}</math> rather than displacements. The  constitutive laws are typically formulated in terms of strain-rate 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  thermo-elasto-visco-plastic model (generally adopted for metals in the  Lagrangian formulation), a thermo-rigid-visco-plastic 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 <math display="inline">\mathbf{\dot{\varepsilon }}=\nabla ^{s}\mathbf{v}</math>. The split  into volumetric and deviatoric parts results in:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\dot{\varepsilon }}=\frac{1}{3}\dot{e}_{vol}\mathbf{I}+\mathbf{\dot{e}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

where <math display="inline">\dot{e}_{vol}</math> and <math display="inline">\mathbf{\dot{e}}</math> 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 elasto-plastic (deviatoric) deformation:

<span id="eq-39"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\dot{e}_{vol}=\mathbf{\nabla }\cdot \mathbf{v}\cong 0    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

so that <math display="inline">\mathbf{\dot{e}}=\mathbf{\dot{\varepsilon }}=\nabla ^{s}\mathbf{v}</math>.

According to the split of the stress tensor introduced in Eq. ([[#eq-28|28]]), different rate-dependent 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, <math display="inline">\mathbf{\dot{e}}</math>. This reduces into the definition of visco-plastic  constitutive models of the form:

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{s}=\eta \mathbf{\,\dot{e}}=2\overline{   \mu   }\mathbf{\dot{e}}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

where <math display="inline">\overline{   \mu   }</math> is the effective viscosity coefficient of the material. The following  additive decomposition of the deviatoric strain rate is assumed:

<span id="eq-41"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\dot{e}}=\mathbf{\dot{e}}^{e}+\mathbf{\dot{e}}^{vp}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

where <math display="inline">\mathbf{\dot{e}}^{e}</math> and <math display="inline">\mathbf{\dot{e}}^{vp}</math> are the elastic and  visco-plastic parts, respectively. In local level FSW, the elastic part of  the strain tensor, <math display="inline">\mathbf{\dot{e}}^{e}</math>, is negligible compared with the  visco-plastic component, <math display="inline">\mathbf{\dot{e}}^{vp}</math>, so that:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{\dot{e}}\equiv \mathbf{\dot{e}}^{vp}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

This means that the total deformation is visco-plastic and a  rigid-visco-plastic model is recovered. The rigid-visco-plastic models can  be interpreted as non-Newtonian laws where the viscosity <math display="inline">\overline{   \mu   }</math> is a nonlinear function of the strain-rate <math display="inline">\mathbf{\dot{e}}.</math> Typical  constitutive equations used for non-Newtonian modeling are listed in Table  constitutive lists. The coefficients <math display="inline">   \mu   _{0}</math> and <math display="inline">   \mu   _{\infty }</math> in these models represent the asymptotic values of the viscosity <math display="inline">\overline{   \mu   }</math> at zero and infinite strain-rates, respectively. In the Newtonian case <math display="inline"> \overline{   \mu   }=   \mu   _{0}=   \mu   _{\infty }=   \mu   .</math>

In the present paper, depending on the definition for <math>\overline{   \mu   }</math>, three models are introduced from the groups (b), (c) and (f):  Norton-Hoff, Sheppard-Wright and Carreau constitutive models.

'''Norton-Hoff model'''     The Norton-Hoff model was originally introduced by Norton <span id='citeF-90'></span>[[#cite-90|[90]]] in  order to describe the unidimensional creep of steel at high temperature, and  extended by Hoff <span id='citeF-70'></span>[[#cite-70|[70]]] to multi-dimensional 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 Norton-Hoff 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, <math> \overset{\cdot }{\bar{\varepsilon }}=\sqrt{2/3}\left\Vert \mathbf{\dot{e}}^{vp}\right\Vert =\sqrt{2/3}\left(\mathbf{\dot{e}}^{vp}:\mathbf{\dot{e}}^{vp}\right)^{1/2}</math>, as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \bar{   \mu }=   \mu   \left(\sqrt{3}\overset{\cdot }{\bar{\varepsilon }}\right)^{(m-1)}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

where <math>0\leq m\left(T\right)\leq 1</math> and <math>   \mu   \left(T\right)</math> are the (temperature dependent) rate sensitivity and  viscosity parameters, respectively (Figure [[#img-17|17]]). Note that:

* <math display="inline">m=0</math> leads to the expression of a rigid plastic constitutive law.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{s}=\sqrt{2}\mu   \mathbf{n}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

where <math display="inline">\mathbf{n=}\dfrac{\mathbf{\dot{\varepsilon }}}{\left\Vert \mathbf{\dot{ \varepsilon }}\right\Vert }=\dfrac{\mathbf{s}}{\left\Vert \mathbf{s}\right\Vert }</math> defines the plastic-flow direction.
* <math display="inline">m=1</math> leads to the linear behavior of a Newtonian fluid.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{s}=2   \mu   \mathbf{\dot{\varepsilon }}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

that is, a linear relationship between stresses and strain-rates where the  constant of proportionality is the viscosity parameter.  

The great advantage of the Norton-Hoff model is its simplicity as it only  depends on two (temperature dependent) parameters.    For the majority of metals, the rate-sensitivity parameter is usually within  the range <math>0.1\leq m\leq 0.2</math> which results in a very non-linear  (non-Newtonian) material behavior (e.g. aluminium and steel alloys).    The nonlinearity of the Norton-Hoff model increases as <math>m</math> diminishes. To  deal with this nonlinearity, <math>m</math> can be computed in a series of steps  starting from the value <math>m=1</math> in the first step and reaching its real value  at step <math>N</math>. This implies that at each time-step <math>i</math> <math><N</math> the problem is  solved with constitutive rate sensitivity <math>m_{i}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  m_{i}=1-\left(1-m\right)\ast \left(\left(1-i\right)/\left(1-N\right) \right)^{m}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

<div id='img-17'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-n-h.png|600px|Stress vs. strain rate affected  by rate sensitivity variation.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 17:''' Stress vs. strain rate affected  by rate sensitivity variation.
|}

'''Sheppard-Wright model'''     The equation used to describe the flow strength of alloys was introduced  initially by Sellars and Tegart <span id='citeF-105'></span>[[#cite-105|[105]]]. Sheppard and Wright <span id='citeF-106'></span>[[#cite-106|[106]]] introduced an alternative rigid visco-plastic formulation  extending this formulation.    According to the Sheppard-Wright model, the effective viscosity coefficient  is defined in terms of effective flow stress <math>\sigma _{e}</math> (i.e. norm of  deviatoric stress) and effective strain rate <math>\overset{\cdot }{\bar{ \varepsilon }}</math> (i.e. norm of deviatoric strain rate) as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \bar{\mu }=\frac{\sigma _{e}}{3\overset{\cdot }{\bar{\varepsilon }}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

where the effective stress <math>\sigma _{e}</math> depends on the strain rate and the  temperature field (below the metals solidus temperature) and it is defined  as

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \sigma _{e}\left(\overset{\cdot }{\bar{\varepsilon }},T\right)=\frac{1}{ \alpha }\sinh ^{-1}\left(\left(\frac{Z}{A}\right)^{1/n}\right)=\frac{1}{ \alpha }\ln \left[\left(\frac{Z}{A}\right)^{1/n}+\sqrt{1+\left(\frac{Z}{A }\right)^{2/n}}\right]   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

<div id='img-18'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-shep.png|600px|Effect of activation energy on  stress vs. strain rate.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 18:''' Effect of activation energy on  stress vs. strain rate.
|}
where <math>Z</math>, the Zener&#8211;Hollomon parameter, represents the  temperature compensated effective strain rate <span id='citeF-121'></span>[[#cite-121|[121]]]

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  Z=\overset{\cdot }{\bar{\varepsilon }}\exp (\frac{Q}{RT})  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

<math>R</math> is the universal gas constant and <math>T</math> is the absolute temperature (in  Kelvin). The material constants <math>\alpha ,</math> <math>A,</math> <math>n</math> and the apparent  activation energy, <math>Q,</math> are derived by fitting experimental response of the  material, and are all independent of temperature (Figure [[#img-18|18]]).    The advantage of the Sheppard-Wright model in comparison with the  Norton-Hoff 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 (<span id='citeF-26'></span>[[#cite-26|[26]]], <span id='citeF-102'></span>[[#cite-102|[102]]]), the effective viscosity coefficient is defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \bar{\mu }=\left[1+\left(\left(\frac{\sigma _{0}}{3\dot{\varepsilon }_{0}\mu _{0}}\right)^{\frac{n}{1-n}}\frac{\sqrt{\frac{2}{3}}\left\Vert   \mathbf{\dot{e}}\right\Vert }{\dot{\varepsilon }_{0}}\right)^{2}\right]^{ \frac{1-n}{2n}}\left(\mu _{0}-\mu _{\infty }\right)+\mu _{\infty }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

<div id='img-19'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-carreau.png|600px|Stress-Strain rate behavior of a  Carreau model with temperature change.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 19:''' Stress-Strain rate behavior of a  Carreau model with temperature change.
|}
where <math>\mu _{0}</math>, <math>\mu  _{\infty }</math>, <math>\dot{\varepsilon }_{0}</math> and <math>n</math> are the saturation viscosities  at zero and infinite limits, the plastic strain rate corresponding to the  yield limit <math>\sigma _{0}</math> and <math>n</math> the exponent, respectively. The yield  limit <math>\sigma _{0}</math> is defined by a Johnson&#8211;Cook power of the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \sigma _{0}=\left\{  \begin{array}{cc}  \sigma _{0,R}\left[1-\left(\frac{T-T_{R}}{T_{M}-T_{R}}\right)^{m}\right] & T<T_{M} \\   0 & T\geq T_{M} \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

where <math>T_{M}</math>, <math>T_{R}</math>, <math>\sigma _{0,R}</math> and <math>m</math> are melting temperature,  room temperature, yield limit at room temperature and Johnson&#8211;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 non-Newtonian constitutive models and thus use  of numerical tricks such as definition of cut-off values on the strain rates.    Figure [[#img-19|19]] shows the stress-strain rate behavior of a Carreau  model for different temperatures.      Neither Carreau and Sheppard-Wright nor Norton Hoff model consider  thermo-elastic 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.

===3.2 Weak form of the mechanical problem===

The mixed <math>\mathbf{u}/p</math> 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 <math>\mathbf{u}/p</math> formulation solves the balance of momentum equation  ([[#eq-30|30]]) and the volumetric part of the constitutive Eq. ( [[#eq-31|31]]) in the weak form. Taking into account that <math>\mathbf{s}\left(  \mathbf{\mathbf{u}}\right)</math> is defined as a function of displacement, the  weak form of the mixed <math>\mathbf{u}/p</math> mechanical problem at global level is  defined over the integration domain <math>\Omega </math> and its boundary <math>\partial  \Omega </math> as:

<span id="eq-52"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{lcc}  \int \limits _{\Omega }\left(\mathbf{s}\left(\mathbf{u}\right):\nabla  ^{s}\delta \mathbf{\mathbf{u}}\right)\mathbf{ }dV+\int \limits _{\Omega  }\left(p\mathbf{\nabla }\cdot \delta \mathbf{\mathbf{u}}\right) dV=W_{mech} &  & \forall \delta \mathbf{v} \\   \int \limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{u}\right)  \mathbf{ }\delta p\right] dV-\int \limits _{\Omega }\left[\left(e^{\theta  }+\dfrac{p}{K}\right)\mathbf{ }\delta p\right] dV=0 &  & \forall  \delta p \end{array}\right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

where <math>\delta \mathbf{\mathbf{u}}</math> and <math>\delta p</math> are the test functions of  the displacement and pressure fields, respectively, while the mechanical  work is defined as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  W_{mech}=\int \limits _{\Omega }\left(\rho _{o}\mathbf{b}\cdot \delta \mathbf{ u}\right)\mathbf{ }dV+\int \limits _{\partial \Omega _{\sigma }}\left(  \mathbf{\bar{t}}\cdot \delta \mathbf{\mathbf{u}}\right)\mathbf{ }dS  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

Definition of the problem in the mixed <math>\mathbf{u/}p</math> format is not  obligatory, however it is convenient for including both the compressible and  incompressible cases. If <math>e^{\theta }\rightarrow 0</math> and the bulk modulus <math> K\rightarrow \infty ,</math> such as in the liquid-like phase, the  incompressibility condition is recovered: <math>\mathbf{\nabla }\cdot \mathbf{u}=0 </math>. In this case, the second Eq. ([[#eq-52|52]]) enforces (in weak  sense) the incompressibility condition as a purely kinematic constraint.    Similarly, according to Eqs. ([[#eq-30|30]]) and ([[#eq-39|39]]),  and taking into account that <math>\mathbf{s}\left(\mathbf{\mathbf{v}}\right)</math>  is defined as a function of the velocity field, the weak form of the mixed <math> \mathbf{v}/p</math> mechanical problem at local level is defined over the  integration domain <math>\Omega </math> and its boundary <math>\partial \Omega </math> as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{lcc}  \int \limits _{\Omega }\left(\mathbf{s}\left(\mathbf{v}\right):\nabla  ^{s}\delta \mathbf{v}\right)\mathbf{ }dV+\int \limits _{\Omega }\left(p \mathbf{\nabla }\cdot \delta \mathbf{\mathbf{v}}\right) dV=W_{mech} &  &  \forall \delta \mathbf{v} \\   \int \limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{v}\right)  \mathbf{ }\delta p\right] dV=0 &  & \forall \delta p \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

where <math>\delta \mathbf{\mathbf{v}}</math> is the test function of the velocity  field.

====3.2.1 Boundary conditions====

For the global level analysis, the boundary <math>\partial \Omega </math> is split into <math>\partial \Omega _{\sigma }</math> and <math>\partial \Omega _{u}</math>, being <math>\partial  \Omega =</math> <math>\partial \Omega _{\sigma }</math> <math>\cup </math> <math>\partial \Omega _{u}</math> such  that the tractions are prescribed on <math>\partial \Omega _{\sigma }</math> while  displacements are specified on <math>\partial \Omega _{u}</math>. The corresponding  Neumann and Dirichlet boundary conditions on <math>\partial \Omega _{\sigma }</math>  and <math>\partial \Omega _{u}</math> are defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{ccc}  \mathbf{\sigma }\cdot \mathbf{n}=\mathbf{\bar{t}} & \hbox{ on} & \partial  \Omega _{\sigma } \\   \mathbf{u}=\mathbf{\bar{u}} & \hbox{ on} & \partial \Omega _{u} \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

where <math>\mathbf{\bar{t}}</math> and <math>\mathbf{\bar{u}}</math> are the prescribed traction  and velocity field, respectively and <math>\mathbf{n}</math> is the unit outward normal  to the surface <math>\partial \Omega _{\sigma }</math>.    Similarly, for the local level analysis, the boundary <math>\partial \Omega </math> is  split into <math>\partial \Omega _{\sigma }</math> and <math>\partial \Omega _{v}</math>, being <math> \partial \Omega =</math> <math>\partial \Omega _{\sigma }</math> <math>\cup </math> <math>\partial \Omega  _{v} </math> such that the velocities are specified on <math>\partial \Omega _{v}</math> and  the corresponding Dirichlet and Neumann boundary conditions on <math>\partial  \Omega _{\sigma }</math> and <math>\partial \Omega _{v}</math> are defined.

===3.3 Discrete weak form of the mechanical problem===

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

* Global form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{lc}  \int \limits _{\Omega }\left(\mathbf{s}_{h}\left(\mathbf{\mathbf{u}}_{h}\right):\nabla ^{s}\delta \mathbf{\mathbf{u}}_{h}\right)\mathbf{ }dV+\int \limits _{\Omega }\left(p_{h} \mathbf{\nabla }\cdot \delta \mathbf{ \mathbf{u}}_{h}\right) dV=W_{mech}^{ext}\left(\delta \mathbf{\mathbf{u}}_{h}\right) & \forall \delta \mathbf{\mathbf{u}}_{h} \\   \int \limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{u}_{h}\right)\mathbf{ }\delta p_{h}\right] dV-\int \limits _{\Omega }\left[ \left(e^{\theta }\left(T_{h}\right)+\dfrac{p_{h}}{K}\right)\mathbf{ }\delta p_{h}\right] dV=0 & \forall \delta p_{h} \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

* Local form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{lc}  \int \limits _{\Omega }\left(\mathbf{s}_{h}\left(\mathbf{\mathbf{v}}_{h}\right):\nabla ^{s}\delta \mathbf{\mathbf{v}}_{h}\right)\mathbf{ }dV+\int \limits _{\Omega }\left(p_{h} \mathbf{\nabla }\cdot \delta \mathbf{ \mathbf{v}}_{h}\right) dV=W_{mech}^{ext}\left(\delta \mathbf{\mathbf{v}}_{h}\right) & \forall \delta \mathbf{\mathbf{v}}_{h} \\   \int \limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{v}_{h}\right)\mathbf{ }\delta p_{h}\right] dV=0 & \forall \delta p_{h} \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

where <math>\mathbf{\mathbf{u}}_{h},</math> <math>\mathbf{\mathbf{v}}_{h},</math> <math>p_{h},</math> <math> \mathbf{s}_{h}</math> and <math>T_{h}</math> are the discrete counterparts of the respective  fields.

==4 Stabilization methods==

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 inf-sup  condition (see first subsection), and avoiding the oscillations arising in  convective dominated problems (see second subsection).

===4.1 Pressure stabilization===

It is known that the incompressibility constraint or under isochoric strain  conditions, the use of an irreducible velocity-based (or displacement-based)  formulation leads to pressure locking situations. It must be pointed out  that when using a J2-plasticity model (in solid mechanics) or one of the  non-Newtonian 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 <math>\mathbf{v}/p</math> formulation is a well-known 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 Ladyzhenskaya-Babuska&#8211;Brezzi (LBB) stability condition <span id='citeF-18'></span>[[#cite-18|[18]]].    The LBB compatibility condition restricts the choice of finite element  spaces <math>V_{h}</math> and <math>P_{h}</math> for velocities (or displacements) and pressure,  respectively. It states that in order to obtain a stable formulation, <math>V_{h}</math>  and <math>P_{h}</math> must be chosen such that the following inf-sup condition is  satisfied:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \underset{\delta p_{h}\in P_{h}}{\inf }\underset{\delta \mathbf{\mathbf{v}}_{h}\in V_{h}}{\sup }\frac{\left(\delta p_{h},\nabla \cdot \delta \mathbf{ \mathbf{v}}_{h}\right)}{\left\Vert \delta p_{h}\right\Vert _{0}\left\Vert  \delta \mathbf{\mathbf{v}}_{h}\right\Vert _{1}}\geq \alpha >0,  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

where parameter <math>\alpha </math> is independent of the mesh size <math>h</math>. If the LBB  compatibility condition is satisfied, there exists a unique solution <math> \mathbf{\mathbf{v}}_{h}\in V_{h}</math> and a <math>p_{h}\in P_{h}</math> (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 <span id='citeF-71'></span>[[#cite-71|[71]]]) or by using  pressure stabilization methods.    Using LBB stable elements, such as the Q2/Q1 (Taylor-Hood 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'' residual-based''. A stabilization method is ''consistent'' if, on convergence, the  solution of the modified problem is exactly the solution of the original  problem.    The sub-grid method was originally introduced in <span id='citeF-73'></span>[[#cite-73|[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 sub-grid (''finer'')  part. While the first one can be resolved by the finite element method, the  latter one cannot. The particular approximation used for the sub-grid 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  sub-grid 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], <span id='citeF-2'></span>[[#cite-2|[2]]], <span id='citeF-37'></span>[[#cite-37|[37]]], [Cervera et
al. 2003]).    To this end, the solution space, <math>\Psi </math>, is split into two parts:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \Psi =\Psi _{h}\oplus \tilde{\Psi }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

<math>\Psi _{h}</math> denotes the finite element solution while <math>\tilde{\Psi }</math>, is the  one that completes <math>\Psi _{h}</math> in <math>\Psi </math> and it is called the '' sub-grid scale'' (<span id='citeF-73'></span>[[#cite-73|[73]]] and <span id='citeF-72'></span>[[#cite-72|[72]]]). Accordingly, the  velocity and pressure fields of the mixed problem are expressed as (Figure  [[#img-20|20]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{v} =\mathbf{v}_{h}+\widetilde{\mathbf{v}} </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|-
| style="text-align: center;" | <math> p =p_{h}   </math>
|}
|}

<div id='img-20'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-sub-scale.png|600px|Finite element and sub-grid scale  parts of the solution U=vpmath]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 20:''' Finite element and sub-grid scale  parts of the solution U=vpmath
|}
where <math>\mathbf{v}_{h}\in </math> <math>\Psi _{h}</math> is the velocity field of the finite element scale and <math> \widetilde{\mathbf{v}}\in \tilde{\Psi }</math> is the enhancement of the velocity  field belonging to the sub-grid scale. Note that no sub-grid scale  contribution is considered for the pressure field. It is reasonable to  assume that the sub-grid velocities <math>\widetilde{\mathbf{v}}</math> are  sufficiently ''small'' compared to <math>\mathbf{v}_{h}</math>. They can be viewed  as a ''high frequency'' perturbation of the finite element field, which  cannot be resolved in <math>\Psi _{h}</math>. It can also be assumed that <math>\widetilde{ \mathbf{v}}</math> and its variation <math>\delta \widetilde{\mathbf{v}}</math> vanish on the  boundary <math>\partial \Omega </math>.    The stress deviator can be similarly decomposed into two different parts:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{s}\left(\mathbf{v}\right)\cong \mathbf{s}_{h}\left(\mathbf{v}_{h}\right)+\widetilde{\mathbf{s}}\left(\widetilde{\mathbf{v}}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

Notice that the stress <math>\mathbf{s}</math> is not a linear function of <math>\mathbf{v}</math>  and, therefore, taking <math>\mathbf{s}_{h}=\mathbf{s}_{h}\left(\mathbf{v}_{h}\right)</math> and <math>\widetilde{\mathbf{s}}=\widetilde{\mathbf{s}}\left(  \widetilde{\mathbf{v}}\right)</math> is a further approximation.    Different choices are available for the approximation of the velocity  sub-grid scale. Two possible options correspond to the ASGS and the OSGS  methods.    Using ASGS method, <math>\widetilde{\mathbf{v}}</math> is approximated as <math>\widetilde{ \mathbf{v}}=</math> <math>\tau _{v}r_{h}</math>, where <math>r_{h}</math> is the residual of the  momentum equation in finite element space.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  r_{h}=r_{h}\left(\mathbf{\mathbf{v}}_{h},p_{h}\right)=\mathbf{\nabla }\cdot \mathbf{s}_{h}+\mathbf{\nabla }p_{h}+\rho _{o}\mathbf{b\simeq \nabla }p_{h}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

where <math>\mathbf{\nabla }\cdot \mathbf{s}_{h}</math> vanishes when using linear  triangular or tetrahedral elements and the body forces, <math>\rho _{o}\mathbf{b}</math> , are negligible. This given:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \widetilde{\mathbf{v}}=\tau _{v}\mathbf{\nabla }p_{h}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 :

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{ccc}  \int\limits _{\Omega }\left(\mathbf{s}_{h}:\nabla ^{s}\delta \mathbf{ \mathbf{v}}_{h}\right)\mathbf{ }dV+\mathbf{\int\limits _{\Omega }}\left( p_{h} \mathbf{\nabla }\cdot \delta \mathbf{\mathbf{v}}_{h}\right) dV=W_{mech}^{ext}\left(\delta \mathbf{v}_{h}\right)&  & \forall \delta   \mathbf{v}_{h} \\   \int\limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{v}_{h}\right)\mathbf{ }\delta p_{h}\right] dV+\int\limits _{\Omega }\tau  _{v}\left(\mathbf{\nabla }\delta p_{h}\cdot \mathbf{\nabla }p_{h}\right) dV=0 &  & \forall \delta p_{h} \end{array}\right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

where <math>\int\limits _{\Omega }\tau _{v}\left(\mathbf{\nabla }\delta  p_{h}\cdot \mathbf{\nabla }p_{h}\right)dV</math> is the resulting stabilization  term to be considered.    A further approximation for the unknown sub-grid space <math>\tilde{\Psi },</math> is  the space orthogonal to the finite element space, referred to hereafter as <math> \Psi _{h}^{\perp }</math>. This means approximating the solution space as <math>\Psi  \simeq \Psi _{h}\oplus \Psi _{h}^{\perp }</math>. The subsequent stabilization  method is called OSGS method, and it has already been successfully applied  to several problems in fluid and solid mechanics. According to OSGS, <math> \widetilde{\mathbf{v}}</math> is approximated as <math>\widetilde{\mathbf{v}}=</math> <math>\tau  _{v}P_{h}^{\perp }\left(r_{h}\right)</math>, where <math>P_{h}^{\perp }</math> is the  orthogonal projection onto the space orthogonal to the finite element space  (Figure [[#img-21|21]]).<div id='img-21'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-osgs.png|600px|Representation of OSGS solution space]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 21:''' Representation of OSGS solution space
|}
The orthogonal projection of a variable <math>\left(\cdot \right)</math> can be  computed as: <math>P_{h}^{\perp }\left(\cdot \right)=\left(\cdot \right) -P_{h}\left(\cdot \right)</math>. The resulting discrete OSGS stabilized weak  form of the mechanical problem is:

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{ccc}  \int\limits _{\Omega }\left(\mathbf{s}_{h}:\nabla ^{s}\delta \mathbf{ \mathbf{v}}_{h}\right)\mathbf{ }dV+\int\limits _{\Omega }\left(p_{h}   \mathbf{\nabla }\cdot \delta \mathbf{\mathbf{v}}_{h}\right) dV=W_{mech}^{ext}\left(\delta \mathbf{v}_{h}\right)&  & \forall \delta   \mathbf{v}_{h} \\   \int\limits _{\Omega }\left[\left(\mathbf{\nabla }\cdot \mathbf{v}_{h}\right)\mathbf{ }\delta p_{h}\right] dV+\int\limits _{\Omega }\tau  _{v}\left[\mathbf{\nabla }\delta p_{h}\cdot \left(\mathbf{\nabla }p_{h}-\Pi _{h}\right)\right]dV=0 &  & \forall \delta p_{h} \\   \int\limits _{\Omega }\left(\mathbf{\nabla }p_{h}\cdot \delta \Pi  _{h}\right)dV-\int\limits _{\Omega }\left(\Pi _{h}\cdot \delta \Pi  _{h}\right)dV=0 &  & \forall \delta \Pi _{h} \end{array}\right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}

where <math>\Pi _{h}=P_{h}\left(\mathbf{\nabla }p_{h}\right)</math> and <math>\delta \Pi  _{h}</math> 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 ([[#eq-65|65]]) is computed  as a function of the orthogonal projection of the residual of the momentum  balance equation as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \int\limits _{\Omega }\tau _{v}\left[\mathbf{\nabla }\delta p_{h}\cdot  P_{h}^{\perp }\left(\mathbf{\nabla }p_{h}\right)\right]dV  </math>
|}
|}

The resulting system is an accurate, stable and consistent formulation for  the solution of the mechanical problem subjected to the incompressibility  constraint. <div id='img-22'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-press-st.png|600px|Pressure field  obtained using ASGS (top) and OSGS (bottom) stabilization method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 22:''' Pressure field  obtained using ASGS (top) and OSGS (bottom) stabilization method.
|}
<div id='img-23'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-conv-1.png|600px|Streamlines obtained using ASGS (top) and OSGS (bottom)  stabilization method.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 23:''' Streamlines obtained using ASGS (top) and OSGS (bottom)  stabilization method.
|}
The stabilization parameter <math>\tau _{v}</math>, is computed as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \tau _{v}=c_{u}\frac{h^{2}}{2\bar{\mu }}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

where <math>h</math> is the element size, <math>c_{u}</math> is a constant and <math>\bar{\mu }</math> is the  parameter defined according to the constitutive model (either effective  viscosity or shear modulus). <math>\bar{\mu }\left(T\right)</math> is usually  temperature dependent and therefore the stabilization term <math>\tau _{v}\left( T\right)</math> is a temperature dependent parameter.    Figure [[#img-22|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 <span id='citeF-6'></span>[[#cite-6|[6]]].    Figure [[#img-23|23]] shows the streamlines. The results  obtained using OSGS are totally unaffected by the boundary conditions.

===4.2 Convective stabilization===

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 well-known 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 <math>10^{1}</math> to <math>10^{3}</math>.  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 sub-grid scale method, the solution space <math>\Theta </math> of the  temperature field is split into two parts: the finite element space and the  space for the sub-grid scale.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \Theta =\Theta _{h}+\tilde{\Theta }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

Therefore, the temperature field can be approximated as

<span id="eq-68"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  T=T_{h}+\tilde{T}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

where <math>T_{h}\in \Theta _{h}</math> is the temperature component of the (coarse)  finite element scale and <math>\tilde{T}\in \tilde{\Theta }</math> is the enhancement of  the temperature field corresponding to the (finer) sub-grid scale. Let us  also consider the corresponding variations <math>\delta T_{h}</math> <math>\in \Theta _{h}</math>  and <math>\delta \tilde{T}</math> <math>\in \tilde{\Theta }</math>, respectively.    Also in this case, we consider two possible choices for the approximation of  the temperature sub-grid scales corresponding to the ASGS and the OSGS  methods.    In the ASGS method, <math>\tilde{T}</math> is approximated as <math>\tilde{T}\mathbf{=}</math> <math> \tau _{\theta }r_{h}^{\theta }</math>, where <math>r_{h}^{\theta }</math> is the residual of  the energy balance equation in the finite element scale.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  r_{h}^{\theta }=\rho _{o}c\left(\frac{\partial T_{h}}{\partial t}+\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) -\mathbf{\nabla \cdot }\left(k \mathbf{\nabla }T_{h}\right)\simeq \rho _{o}c\left(\frac{\partial T_{h}}{ \partial t}+\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}

where the diffusive term, <math>\mathbf{\nabla \cdot }\left(k\mathbf{\nabla }T_{h}\right),</math> 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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \int\limits _{\Omega }\tau _{\theta }\left[\rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }\left(\delta T_{h}\right)\right)\right] \left[ \rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) \right] dV  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

to the Galerkin weak form of the problem. This given:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{ccc}  \int\limits _{\Omega }\left[\rho _{o}c\left(\dfrac{\partial T_{h}}{ \partial t}+\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) \delta T_{h} \right] dV+\int\limits _{\Omega } \left(k\nabla T_{h}\right)\mathbf{ \cdot }\nabla \left(\delta T_{h}\right)dV  &  &  \\   +\int\limits _{\Omega }\tau _{\theta }\left[\rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }\left(\delta T_{h}\right)\right)\right] \left[ \rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) \right] dV=W_{ther}^{ext}\left(\delta T_{h}\right)&  & \forall \delta T_{h} \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 <math>\Theta  _{h}^{\perp }</math> can be adopted for the approximation of the unknown sub-grid  space <math>\tilde{\Theta }</math>. This means that the solution space is approximated  as <math>\Theta \simeq \Theta _{h}\oplus \Theta _{h}^{\perp }</math>. According to the  OSGS method, <math>\tilde{T}</math> is defined as <math>\tilde{T}\mathbf{=}</math> <math>\tau _{\theta  }P_{h}^{\perp }\left(r_{h}^{\theta }\right)</math>.    Taking into account that <math>P_{h}^{\perp }\left(\cdot \right)=\left(\cdot  \right)-P_{h}\left(\cdot \right)</math>, the resulting discrete OSGS stabilized  weak form of the thermal problem for linear elements simplifies to:

<span id="eq-72"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{ccc}  \int\limits _{\Omega }\left[\rho _{o}c\left(\dfrac{\partial T_{h}}{ \partial t}+\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right) \delta T_{h} \right] dV+\int\limits _{\Omega }\left(k\nabla T_{h}\right)\mathbf{\cdot  \nabla }\left(\delta T_{h}\right)dV  &  &  \\   +\int\limits _{\Omega }\tau _{\theta }\left[\rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }\left(\delta T_{h}\right)\right)\right] \left[ \rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right)-\Pi  _{h}^{\theta } \right]dV=W_{ther}^{ext}\left(\delta T_{h}\right)&  &   \forall \delta T_{h} \\   \int\limits _{\Omega }\left[\left(\mathbf{v}_{h}\cdot \mathbf{\nabla }T_{h}\right)\hbox{ }\delta \Pi _{h}^{\theta } \right]dV-\int\limits _{ \Omega }\left[\Pi _{h}^{\theta }\cdot \hbox{ }\delta \Pi _{h}^{\theta }  \right]dV=0 &  & \forall \delta \Pi _{h}^{\theta } \end{array}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
|}

where <math>\Pi _{h}^{\theta }</math> is the smooth projection of the convective term.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \Pi _{h}^{\theta }=P_{h}\left[\rho _{o}c\left(\mathbf{v}_{h}\cdot \mathbf{ \nabla }T_{h}\right) \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
|}

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

<span id="eq-74"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \tau _{\theta }=\left(c_{1}\frac{k}{h^{2}}+c_{2}\frac{\left\Vert \mathbf{v}_{h}\right\Vert }{h}\right)^{-1}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
|}

where <math>h</math> is the element size, <math>c_{1}</math> and <math>c_{2}</math> are algorithmic  constants. <div id='img-24'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-temp.png|600px|Comparison between  the temperature distribution at different times (10, 40 and 70) and at the  steady state.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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 [[#img-24|24]] in a parallel flow passing  through a tube.

==5 Apropos kinematic framework==

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, <math>\mathbf{X}</math>, coincides with the mesh  nodes. The material and local time derivative of the temperature concide:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \dfrac{dT\left(\mathbf{X},t\right)}{dt}=\dfrac{\partial T\left(\mathbf{X},t\right)}{\partial t}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \dfrac{dT\left(\mathbf{x},t\right)}{dt}=\dfrac{\partial T\left(\mathbf{x},t\right)}{\partial t}+\dfrac{\partial T\left(\mathbf{x},t\right)}{ \partial \mathbf{x}}\dfrac{d\mathbf{x}}{dt}=\dfrac{\partial T\left(\mathbf{x },t\right)}{\partial t}+\mathbf{v}\left(\mathbf{x},t\right)\cdot \nabla _{ \mathbf{x}}T\left(\mathbf{x},t\right) </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
|}

where <math>\mathbf{v}\left(\mathbf{x},t\right)</math> is the spatial velocity and <math> \nabla _{\mathbf{x}}T\left(\mathbf{x},t\right)</math> 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 <math>\mathbf{c}</math>. It is  defined as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}  \dfrac{dT\left(\mathbf{\chi },t\right)}{dt} &\cong &\dfrac{\partial  T\left(\mathbf{\chi },t\right)}{\partial t}+\dfrac{\partial T\left(  \mathbf{\chi },t\right)}{\partial \mathbf{\chi }}\dfrac{d\mathbf{\chi }}{dt} =  \\ &&\dfrac{\partial T\left(\mathbf{\chi },t\right)}{\partial t}+\mathbf{c}\cdot \nabla _{\mathbf{\chi }}T\left(\mathbf{\chi },t\right) \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
|}

where the relative velocity is <math>\mathbf{c}=\mathbf{v}\left(\mathbf{\chi },t\right)-\mathbf{v}_{mesh}</math>. <math>\mathbf{v}\left(\mathbf{\chi },t\right)</math>  and <math>\mathbf{v}_{mesh}</math> are the material and the mesh velocities,  respectively .    In this case, the movement of the mesh must be calculated at each time-step.  In general, an ad-hoc methodology is necessary to compute the position of  the mesh at each time-step 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 work-piece and the back-plate.    Taking into account that, during the welding process, the pin is rotating  with a very high speed (e.g. 50-5000 RPM, depending on the work-piece  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 strain-rates. As a consequence,  a fine mesh and a continuous re-meshing 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 cross-section. 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 re-defined at each time-step 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 extent of the material deformation varies in different regions of  the analysis domain.
* Pins used in practice are not necessarily cylindrical (Figure [[#img-5|5]]).<div id='img-25'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-kin-fra.png|600px|The geometry  discretization]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 25:''' The geometry  discretization
|}
* Boundary conditions must be applied conveniently according to the  movement of the pin.
* Re-definition of the integration domain and re-meshing is preferably  to be avoided.  

In the proposed methodology, the original geometry is divided into three  distinct zones (Figure [[#img-25|25]]) which are: the pin, the work-piece and  the stir-zone (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.

* The pin is assumed to be rigid and it is defined in the Lagrangian  framework as at each time-step of the analysis, the mesh moves according to  the rotation of the pin.
* The work-piece, excluding the stir-zone close to the pin, is modeled  in the Eulerian framework, as in this area the domain neither changes its  shape nor contains moving boundaries. The fixed mesh in Eulerian zone  enables the application of boundary conditions in the inflow and outflow of  the work-piece.
* In the stir-zone, the material is subjected to extremely high strains.  Thus, the ALE kinematic framework is chosen to study the stir-zone. The  stir-zone is modeled as a circular region around the pin. The size of this  area strongly depends on the viscosity and thermal diffusivity of the  material. The key idea consists in using the so-called mesh sliding  techniques, meaning that the mesh rigidly rotates at each time-step  according to the pin movement, decoupling the material motion from the  motion of the mesh. Consequently, the integration domain moves in such a way  that the inner boundary of the stir-zone is connected to the contour surface  of the pin. In the ALE approach, the calculation of the mesh velocity can be  a difficulty while in the proposed strategy, this is not an issue. When  studying a FSW process, the mesh velocity can be prescribed according to the  pin rotation as:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{v}_{mesh}=\mathbf{\varpi \times r}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
|}

where <math display="inline">\mathbf{\varpi }</math> is the angular velocity of the pin and <math display="inline">\mathbf{r}</math>  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 <math>\mathbf{v}_{mesh}</math> the mesh  velocity, the convective term can be written as <math>\rho _{o}c\left[\left(  \mathbf{v}-\mathbf{v}_{mesh}\right)\cdot \mathbf{\nabla }T\right]</math>, the  convective term vanishes in the Lagrangian framework ( <math>\mathbf{v=v}_{mesh}</math> ) and simplifies to <math>\rho _{o}c</math> <math>\mathbf{v}\cdot \mathbf{\nabla }T</math> in the  Eulerian zone ( <math>\mathbf{v}_{mesh}=\mathbf{0}</math>). Table [[#table-2|2]]  summarizes the computational framework together with the solution hypotheses  for the pin, the work-piece and the stir-zone.         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 work-piece and the ALE mesh used for the stir-zone must be  considered. On the other hand, the frictional contact between the stir-zone  and the pin (see section 2.5) must be defined.    The coupling at the interface between the ALE (stir-zone) and the Eulerian  (work-piece) parts requires a special treatment as the mesh is fixed on one  side (Eulerian) and it moves on the other side (ALE). Physically, work-piece  and stir-zone are not separated, they belong to the same metal sheet even if  in the numerical model, a relative movement of the two computational  sub-domains exists. Therefore, the objective is to get continuous fields for  all the state variables <math>\mathbf{v,}</math> <math>p</math> and <math>T</math> at the interface between  the two zones. The coupling is performed using a ''node-to-node link'' or   ''node-to-face contact ''approach.<div id='img-26'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-master-slave.png|600px|Node-to-node link approach.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 26:''' Node-to-node link approach.
|}
In the first and simpler case, each node of the mesh on the interface is  duplicated (Figure [[#img-26|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  node-to-node connections at the interface. This leads to an ad-hoc  assembling procedure, where the contributions (elemental residuals and  tangent matrices) of the stir-zone nodes are summed to the corresponding  contributions of the work-piece nodes. This approach requires a coincident  mesh at the interface. The time-step is chosen such that the two meshes (ALE  and Eulerian) are synchronized to keep the Courant number, <math>Cu=1</math> at the  sliding interface.    The node-to-node 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 time-step, the ALE mesh at the interface slides from one node to another. This requires  coincident and equi-spaced mesh at the contact surfaces.    If the surface meshes are not coincident, a node-to-face 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 stir-zone is projected on the work-piece  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, <math>T</math>, <math>p</math>, <math> \mathbf{v}</math>.    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 time-step, the contact elements need to be  regenerated according to a (rather time-consuming) closest-point-projection  algorithm.

==6 Comparison between a Solid Mechanics and a Fluid Mechanics approach  for the numerical simulation of FSW processes at local level==

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 ad-hoc Apropos Kinematic Framework  (AKF) which makes use of an efficient combination of Lagrangian (tool),  Eulerian (work-pieces) and ALE (stirring zone) descriptions for the  different computational sub-domains, within either a Solid Mechanics or a  Fluid Mechanics approach, is the most suitable option (<span id='citeF-22'></span>[[#cite-22|[22]]],   <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-50'></span>[[#cite-50|[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 (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]). The main features of each one of those approaches, pointing out which are  their main advantages and drawbacks, will be presented.

===6.1 Coupled thermo-mechanical problem===

====6.1.1 Solid Mechanics approach====

In the Solid Mechanics approach, the governing equations of the coupled  thermo-mechanical problem are written in terms of the displacement vector  field and temperature scalar field, which are taken as primary variables ( <span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]).    The original coupled thermo-mechanical 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 thermo-mechanical problem leads to a Product Formula Algorithm  (PFA), where the original coupled thermomechanical problem is solved by  means of a staggered time-marching 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.

====6.1.2 Fluid Mechanics approach====

In the Fluid Mechanics approach, the governing equations of the (fully  incompressible) coupled thermo-mechanical 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 (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-50'></span>[[#cite-50|[50]]]).    The original coupled thermo-mechanical 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 thermo-mechanical problem leads to a Product Formula Algorithm  (PFA), where the original coupled thermomechanical problem is solved by  means of a staggered time-marching 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.

===6.2 FEA of the coupled thermomechanical problem===

====6.2.1 Solid Mechanics approach====

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 quasi-incompressibility constraint, the pressure is  considered constant over the element and it is computed only at a central  integration point (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]).

====6.2.2 Fluid Mechanics approach====

In the Fluid Mechanics approach, a sub-grid scale stabilized mixed linear  velocity/pressure/ temperature finite element formulation is used. The  Orthogonal Sub-grid 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 (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-50'></span>[[#cite-50|[50]]], <span id='citeF-6'></span>[[#cite-6|[6]]]).

===6.3 Apropos Kinematic Framework (AKF)===

In order to introduce an AKF suitable for both a Solid Mechanics and a Fluid  Mechanics framework, let us introduce the following computational  sub-domains for the plates. Taking into account the radial distance to the  rotation axis of the tool, three different zones or computational  sub-domains, as shown in Figure [[#img-27|27]], are introduced (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]).    In the closest zone around to the tool (zone 1), the stir-zone, the material  is subjected to extremely high strains. The stir-zone is modeled as a  circular (2D) or cylindrical (3D) region around the tool and it includes the  Thermo-Mechanically Affected Zone (TMAZ). The size of the stir-zone strongly  depends on the viscosity and thermal diffusivity of the material. Due to the  extremely high deformations, the use of a Lagrangian description would lead  very quickly to mesh entanglement. Thus, in this region, an ALE description  is a most suitable option. On the other hand, the use of an ALE description  allows to bypass problems related to moving boundaries which would arise for  non-axisymmetric pins if a Eulerian description is used. The mesh velocity  in the stir-zone is prescribed to the rotational velocity of the tool.<div id='img-27'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-sf1.png|600px|Computational sub-domains of the  plates: Stir-zone (zone 1), transition zone (zone 2), and work-pieces zone  (zone 3)]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 27:''' Computational sub-domains of the  plates: Stir-zone (zone 1), transition zone (zone 2), and work-pieces 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 rigid-thermal body and a Lagrangian description is  used in both the Solid Mechanics and Fluid Mechanics approaches.

====6.3.1 Solid Mechanics approach====

In the Solid Mechanics approach, an ALE description is applied to each one  of the three computational sub-domains defined in the Figure [[#img-27|27]]. In  the stir-zone 1, the mesh velocity is prescribed to the rotational velocity  of the tool. In the work-pieces 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 sub-domain, the rotational tool velocity at the  interface with the zone 1, and zero velocity at the interface with the zone  3 (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[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 <span id='citeF-25'></span>[[#cite-25|[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:

<span id="eq-79"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}  \left.\dfrac{\partial \rho }{\partial t}\right\vert _{\chi }+\mathbf{c}\cdot \nabla \rho +\rho \nabla \cdot \mathbf{v} &=&0   \\ \rho \left(\left.\dfrac{\partial \mathbf{v}}{\partial t}\right\vert _{\chi  }+\mathbf{c}\cdot \nabla \mathbf{v}\right)&=&\nabla \cdot \mathbf{\sigma }+\rho \mathbf{b}  \\ \rho c_{p}\left(\left.\dfrac{\partial T}{\partial t}\right\vert _{\chi }+ \mathbf{c}\cdot \nabla T\right)&=&D_{mech}-\nabla \cdot \mathbf{q}   \end{array}</math>
|}
|}

where <math>\rho </math> is the spatial mass density, <math>c_{p}</math> is the specific heat  capacity, <math>\mathbf{v}</math> is the material velocity vector field, <math>\mathbf{c}</math>  is the convective velocity vector field, defined as the difference between  the material velocity and the mesh velocity vector fields, <math>\mathbf{\sigma }</math>  is the Cauchy stress tensor field, <math>\mathbf{b}</math> is the body forces vector  field per unit of mass, <math>T</math> is the temperature field, <math>D_{mech}</math> is the  mechanical dissipation rate per unit of spatial volumen, <math>\mathbf{q}</math> is the  heat flux vector per unit of spatial surface, and <math>\mathbf{\chi }</math> 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 ([[#eq-79|79]]) is solved using  an operator split procedure, which yields to a product formula algorithm  (PFA) (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-16'></span>[[#cite-16|[16]]], <span id='citeF-17'></span>[[#cite-17|[17]]], <span id='citeF-25'></span>[[#cite-25|[25]]], <span id='citeF-53'></span>[[#cite-53|[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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}  \left(\rho \mathbf{\ddot{u},}\delta \mathbf{u}\right)+\left(\mathbf{ \sigma ,\nabla }\delta \mathbf{u}\right)&=&\left(\rho \mathbf{b,}\delta   \mathbf{u}\right)+\left(\mathbf{t,}\delta \mathbf{u}\right)_{\partial  _{\sigma }\Omega } \\ \left(\rho \mathbf{c}_{p}\dot{T}\mathbf{,}\delta T\right)+\left(k\nabla T \mathbf{,\nabla }\delta T\right)&=&\left(D_{mech}\mathbf{,}\delta T\right) -\left(\bar{q}\mathbf{,}\delta T\right)_{\partial _{q}\Omega }  \\ &&-\left(q_{cond}\mathbf{,}\delta T\right)_{\partial _{cond}\Omega  }-\left(q_{conv}+q_{rad}\mathbf{,}\delta T\right)_{\partial _{env}\Omega }   \end{array}</math>
|}
|}

where <math>\delta \mathbf{u}</math> and <math>\delta T</math> are the test functions for the  displacement vector field and temperature field, respectively, <math>\mathbf{t}</math>  is the prescribed traction vector field on the boundary <math>\partial _{\sigma  }\Omega </math>, <math>\bar{q}</math> is the prescribed normal heat flux on <math>\partial  _{q}\Omega </math>, <math>q_{cond}</math> is the conduction heat fluxes on <math>\partial  _{cond}\Omega </math>, <math>q_{conv}</math> and <math>q_{rad}</math> are the convection and radiation  normal heat fluxes on <math>\partial _{env}\Omega </math>, respectively, and <math>k</math> 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 sub-steps. 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.

====6.3.2 Fluid Mechanics approach====

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  stir-zone 1, where an ALE formulation is used, and the other one is linked  to the work-pieces 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 node-to-node (NTN) direct contact link approach.  This approach requires a coincident and equi-spaced 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 <math>Cu=1</math>. 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 equi-spaced, the meshes cannot be synchronized  and a more general node-to-segment (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  (work-pieces, zone 3) and ALE (stirring zone, stir zone 1) descriptions for  the different computational sub-domains, 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 ( <span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-50'></span>[[#cite-50|[50]]], <span id='citeF-6'></span>[[#cite-6|[6]]]):

<ol>

<li>The flow is assumed to be fully incompressible as a result of  considering a rigid-plastic constitutive model with isochoric plastic  behaviour (J2 plasticity), and considering negligible the volumetric thermal  strains. Then, the spatial divergence of the velocity is zero, the spatial  density is constant, i.e. <math display="inline">\rho =\rho _{0}</math>, and the deviatoric part of the  deformation rate is computed as the symmetric part of the spatial velocity  gradient.    </li>
<li>The problem can be analyzed as quasi-static. The inertial term in the  momentum balance equation is assumed to be negligible due to the very high  viscosity of the material. The material flow is characterized by very low  Reynolds numbers (Re<<1).  </li>

</ol>

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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}  \nabla \cdot \mathbf{v} &=&0 \\ \mathbf{0} &=&\nabla \cdot \mathbf{s}+\nabla p+\rho _{0}\mathbf{b}  \\ \rho _{0}c_{p}\left(\left.\dfrac{\partial T}{\partial t}\right\vert _{\chi  }+\mathbf{c}\cdot \nabla T\right)&=&D_{mech}-\nabla \cdot \mathbf{q}   \end{array}</math>
|}
|}

where <math>\rho _{0}</math> is the reference mass density, <math>c_{p}</math> is the specific  heat capacity, <math>\mathbf{v}</math> is the material velocity vector field, <math>\mathbf{c }</math> is the convective velocity vector field, defined as the difference  between the material velocity and the mesh velocity vector fields, <math>\mathbf{s }</math> is the deviatoric part of the Cauchy stress tensor field, <math>p</math> is the  pressure field, <math>\mathbf{b}</math> is the body forces vector field per unit of  mass, <math>T</math> is the temperature field, <math>D_{mech}</math> is the mechanical dissipation  rate per unit of spatial volumen, <math>\mathbf{q}</math> is the heat flux vector per  unit of spatial surface, and <math>\mathbf{\chi }</math> 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 quasi-static problem, and energy balance equations takes the  form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\begin{array}{l}  \left(\nabla \cdot \mathbf{v,}\delta p\right)&=&0 \\ \left(\mathbf{s,}\nabla \delta \mathbf{v}\right)+\left(p\mathbf{,}\nabla  \cdot \delta \mathbf{v}\right)&=&W_{mech}  \\ \left(\rho _{0}c_{p}\left.\dfrac{\partial T}{\partial t}\right\vert _{\chi  },\delta T\right)+\left(\rho _{0}c_{p}\mathbf{c}\cdot \nabla T,\delta  T\right)+\left(k\nabla T,\nabla \delta T\right)&=&W_{ther}   \end{array}</math>
|}
|}

where

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  W_{mech} =\left(\rho _{0}\mathbf{b,}\delta \mathbf{v}\right)+\left(  \mathbf{t,}\delta \mathbf{v}\right)_{\partial _{\sigma }\Omega } </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
|-
| style="text-align: center;" | <math> W_{ther} =\left(D_{mech}\mathbf{,}\delta T\right)-\left(\bar{q}\mathbf{, }\delta T\right)_{\partial _{q}\Omega }-\left(q_{cond}\mathbf{,}\delta  T\right)_{\partial _{cond}\Omega }-\left(q_{conv}+q_{rad}\mathbf{,}\delta  T\right)_{\partial _{env}\Omega }   </math>
|}
|}

where <math>\delta p</math>, <math>\delta \mathbf{v}</math> and <math>\delta T</math> are the test functions  for the pressure field, displacement vector field and temperature field,  respectively, <math>\mathbf{t}</math> is the prescribed traction vector field on the  boundary <math>\partial _{\sigma }\Omega </math>, <math>\bar{q}</math> is the prescribed normal  heat flux on <math>\partial _{q}\Omega </math>, <math>q_{cond}</math> is the conduction normal  heat flux on <math>\partial _{cond}\Omega </math>, <math>q_{conv}</math> and <math>q_{rad}</math> are the  convection and radiation normal heat fluxes on <math>\partial _{env}\Omega </math>,  respectively, and <math>k</math> is the isotropic thermal conductivity parameter.    Figure [[#img-28|28]] shows the evolution of the temperature computed by the  Solid Mechanics and Fluid Mechanics approaches at three different points   <span id='citeF-22'></span>[[#cite-22|[22]]]. The frequency of the temperatures at points 1 and 2,  both of them rotating with the tool, is controlled by the rotational  velocity of the tool. Both approaches essentially provide the same results  and further investigations would be still needed to understand the small  differences shown by the two approaches. From a computational point of view,  the Fluid Mechanics approach is more efficient than the Solid Mechanics one.  On the other hand, as the Fluid Mechanics approach uses a rigid-plastic  constitutive model, residual stresses cannot be computed using a  straightforward procedure.<div id='img-28'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-sf2.png|600px|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]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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
|}

===6.4 Thermo-mechanical constitutive models===

In both the Solid Mechanics and Fluid Mechanics approaches, the tool is  modeled as a rigid body, using a thermo-rigid constitutive model, and the  plates can be modeled using a thermo-visco-plastic constitutive model, e.g.  a Norton-Hoff 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 Norton-Hoff constitutive model for the plates, the deviatoric  component of the Cauchy stress tensor takes the form:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{s}=2\mu \left(T\right)\mathbf{\dot{e}}\left(\sqrt{3}\sqrt{\frac{2 }{3}\mathbf{\dot{e}:\dot{e}}}\right)^{m\left(T\right)-1}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
|}

where <math>\mathbf{\dot{e}=}\mathrm{dev}\left(\nabla ^{s}\mathbf{v}\right)</math> is  the deviatoric component of the deformation rate, <math>\mu \left(T\right)</math> is  the temperature dependent viscosity parameter, and <math>m\left(T\right)</math> is a  temperature dependent rate sensitivity parameter.    The mechanical dissipation rate per unit of spatial volume is computed as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  D_{mech}=\gamma _{heat}\mathbf{s}:\mathbf{\dot{e}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
|}

where <math>\gamma _{heat}\approx 0.9</math> is a parameter that represents the  fraction of the total plastic dissipation rate converted into heat.

====6.4.1 Solid Mechanics approach====

In the Solid Mechanics approach, two different constitutive models can be  considered (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]). A first option would be to consider a thermo-rigid-plastic model, e.g.  the Norton-Hoff model introduced above, to characterize the deviatoric  component of the Cauchy stress tensor, and a thermo-elastic 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,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \dot{p}=K\left(T\right)\nabla \cdot \mathbf{v}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
|}

where <math>K\left(T\right)</math> is the temperature dependent bulk modulus. Note  that with this constitutive model, residual stresses can not be computed  using a straightforward procedure.    Alternatively, a thermo-elasto-visco-plastic constitutive model suitable for  finite deformations can be used <span id='citeF-94'></span>[[#cite-94|[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.

====6.4.2 Fluid Mechanics approach====

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 two-step 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.

===6.5 Thermo-mechanical contact===

Two different thermo-mechanical contact approaches are considered.

====6.5.1 Solid Mechanics approach====

In the Solid Mechanics approach, full-stick thermo-mechanical contact  conditions are considered between the tool and the work-pieces. Thus, the  displacement and temperature fields are continuous through the contact  interface between the tool and the work-pieces (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]]).

====6.5.2 Fluid Mechanics approach====

In the Fluid Mechanics approach, a thermo-frictional contact model is  considered (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]],   <span id='citeF-50'></span>[[#cite-50|[50]]]). Using a Norton thermo-frictional contact model, the  tangential component of the traction vector at the contact interface between  the tool and the work-pieces is given by,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{t}_{T}=a\left(T\right)\left\Vert \Delta \mathbf{v}_{T}\right\Vert  ^{q-1}\Delta \mathbf{v}_{T}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
|}

where <math>\mathbf{t}_{T}</math> and <math>\Delta \mathbf{v}_{T}</math> 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 work-pieces, <math>a\left(T\right)</math> is the temperature dependent  material consistency parameter, and <math>0\leq q\leq 1</math> is the strain rate  sensitivity parameter.    The heat flux generated by friction at the contact interface between the  tool and the work-pieces is split into a part absorbed by the tool, denoted  as <math>q_{frict}^{tool}</math>, and a part absorbed by the work-pieces, denoted as <math> q_{frict}^{wp}</math>, given, respectively, by (<span id='citeF-22'></span>[[#cite-22|[22]]], <span id='citeF-23'></span>[[#cite-23|[23]]], <span id='citeF-24'></span>[[#cite-24|[24]]], <span id='citeF-50'></span>[[#cite-50|[50]]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{c}  q_{frict}^{tool}=\vartheta _{frict}^{tool}\mathbf{t}_{T}\cdot \Delta \mathbf{ v}_{T}=\vartheta _{frict}^{tool}a\left(T\right)\left\Vert \Delta \mathbf{v}_{T}\right\Vert ^{q+1} \\   q_{frict}^{wp}=\vartheta _{frict}^{wp}\mathbf{t}_{T}\cdot \Delta \mathbf{v}_{T}=\vartheta _{frict}^{wp}a\left(T\right)\left\Vert \Delta \mathbf{v}_{T}\right\Vert ^{q+1} \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
|}

where the amount of heat absorbed by the tool, <math>\vartheta _{frict}^{tool}</math>,  and by the work-pieces, <math>\vartheta _{frict}^{wp}</math>, depends on the thermal  diffusivity of the two materials in contact, and take the form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{c}  \vartheta _{frict}^{tool}=\dfrac{\alpha _{dif}^{tool}}{\alpha  _{dif}^{tool}+\alpha _{dif}^{wp}} \\   \vartheta _{frict}^{wp}=\dfrac{\alpha _{dif}^{wp}}{\alpha  _{dif}^{tool}+\alpha _{dif}^{wp}} \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
|}

where the thermal diffusivity parameter is given by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \alpha _{dif}=\frac{k}{\rho c_{p}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
|}

Alternatively, as a limit case, full-stick thermo-mechanical contact  conditions between the tool and the work-pieces can be also considered. In  this case the velocity and temperature fields are continuous through the  contact interface between the tool and the work-pieces.

===6.6 Concluding remarks===

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 quasi-incompressible  thermo-elasto-plastic material model (a thermo-rigid-plastic 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 thermo-rigid-plastic 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 re-meshing 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 re-meshing 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  thermo-elasto-plastic consitutive model, is used.

==7 Particle tracing technique==

One of the main issues in the study of FSW is heat generation. In the FSW  process, welding is achieved by the heat generated due to friction and the  material mixing/stirring process. The heat generated must be enough to allow  for the material to flow and to obtain a deep heat affected zone. The  visualization of the material flow is very useful to understand its behavior  during the weld (Figure [[#img-29|29]]). A method assessing the  quality of the created weld by visualization of the joint pattern is  advantageous. It can be used to have a pre-knowledge of the appropriate  process parameters.<div id='img-29'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-weld-quality-deep.png|600px|Material flow  visualization in horizontal and vertical directions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 29:''' Material flow  visualization in horizontal and vertical directions.
|}
In this paper, a numerical particle tracing technology is introduced to  study the extent of material stirring during the FSW process and to study  the weld quality.    Using the Lagrangian framework, the solution obtained at the mesh points  represents the material movement: points of the mesh are material particles.  Therefore, no special technique for the tracking of the material is  necessary. A particle tracing technique needs to be used when the solution  of the problem is performed in ALE or Eulerian kinematic frameworks.    Due to the ALE framework of the finite element analysis used, the motion of  the mesh is not necessarily tied to the motion of the material. During the  analysis, a material particle moves through the mesh and at different time,  it is located inside different elements. To observe material movement around  the pin, it is necessary to construct and analyze material particle  trajectories. This is possible with the use of a particle tracing method to  follow the motion of material points. This method can be naturally applied  to the study of the material flow in the welding process.<div id='img-30'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-particle-distribution.png|600px|Material points distribution.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''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  post-process level must be solved (Figure particle distribution). Each  particle position is integrated from:

<span id="eq-91"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \frac{D\left(\mathbf{X}\left(t\right)\right)}{Dt}=\mathbf{V}\left(  \mathbf{X}\left(t\right),t\right)   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
|}

with the initial condition

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{X}\left(t=0\right)=\mathbf{X}_{0}  </math>
|}
|}

<div id='img-31'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-tracer.png|600px|Tracer position in time]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 31:''' Tracer position in time
|}
where <math>\mathbf{X}\left(t\right)</math> is  the position of the material points at time <math>t</math> and <math>\mathbf{V}\left(  \mathbf{X}\left(t\right),t\right)</math> is the velocity of the tracer in the  position <math>\mathbf{X}\left(t\right)</math> and time <math>t</math>. 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 [[#img-31|31]]).    To solve Eq. ([[#eq-91|91]]), there exist a large amount of integration  techniques ranging from the simple first order Backward Euler (BE) scheme to  higher order Runge-Kutta schemes. Among them three well-known methods are  chosen and compared: Backward Euler with Sub-stepping (BES), Fourth order  Runge-Kutta (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 stir-zone leading to a better  understanding of the welding process.

===7.1 BES method===

According to the first-order accurate BES method, the particle position at  time-step <math>n+1</math> is computed in <math>k</math> sub-steps as:

<span id="eq-92"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{c}  \mathbf{X}_{n+1}\cong \mathbf{X}_{n}+\sum \limits _{i=1}^{k}\mathbf{V}\left(  \mathbf{X}_{n+1}^{i},t_{n}+\sum \limits _{j=1}^{i}\delta t^{j}\right)\delta  t^{i}=\mathbf{X}_{n}+\sum \limits _{i=1}^{k}\mathbf{V}_{n+1}^{i}\delta t^{i}  \\   \dfrac{\Delta t}{k}=\delta t^{i} \end{array}\right.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
|}

where <math>\Delta t</math> is the time-step and <math>\delta t^{i}</math> is the <math>i^{th}</math>  sub-step.    The velocity of the tracer at the current sub-step <math>i</math>, <math>\mathbf{V}_{n+1}^{i} </math>, is computed as:

* The tracer position <math display="inline">\mathbf{X}_{n+1}^{i}=\mathbf{X}\left( t_{n}+\sum \limits _{j=1}^{i}\delta t^{j}\right)</math> is updated according to the  previous sub-step;
* The tracer particle is located within a specific FE by a corresponding  searching algorithm;
* The tracer velocity is interpolated from the nodal velocities of the  FE identified according to the sub-stepping time, <math display="inline">\Upsilon _{t}:</math>

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \left\{  \begin{array}{cc}  \mathbf{v}_{n+1}^{i}\left(t\right)=\left(1-\Upsilon _{t}\right)\mathbf{v}_{n}+\Upsilon _{t}\mathbf{v}_{n+1} &  \\   \Upsilon _{t}=\left(\dfrac{t-t_{n}}{t_{n+1}-t_{n}}\right)& t_{n}<t<t_{n+1} \end{array}\right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
|}

Even if the accuracy of the method is enhanced by a sub-stepping technique,  it is still first-order accurate <math>O\left(\Delta t\right).</math>

===7.2 BFECC method===

<div id='img-32'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-bfecc4.png|600px|The steps of the BFECC method for a  rotational velocity field]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 32:''' The steps of the BFECC method for a  rotational velocity field
|}
According to the third-order accurate, <math>O\left(\Delta t^{3}\right)</math>, BFECC  method, the error arising from backward and forward advection of the  particle position is used to modify the initial position.    Denoting with <math>A\left(\mathbf{\cdot }\right)=\left(\mathbf{\cdot }\right) +\int _{t_{n}}^{t_{n+1}}\mathbf{V}\left(\mathbf{\cdot },t\right)\mathrm{d}t, </math> the first order up-winding integration operator, BFECC consists of a  prediction/correction procedure: the initial position of the tracer, <math> \mathbf{X}_{n}=\mathbf{X}\left(t_{n}\right)</math>, is advected forward to <math> \mathbf{\hat{X}}_{n+1}=A\left(\mathbf{X}_{n}\right)</math>, and then backward to <math>\mathbf{\hat{X}}_{n}=A^{R}\left(\mathbf{\hat{X}}_{n+1}\right)</math>, <math> A^{R}\left(\mathbf{\cdot }\right)</math> being the backward advection operator.  The difference between <math>\mathbf{\hat{X}}_{n}</math> and the initial position  defines the error, <math>\mathbf{\epsilon =}\dfrac{\mathbf{X}_{n}-\mathbf{\hat{X}}_{n}}{2}</math>, of the time integration scheme.    This error is used later to correct the starting position, <math>\mathbf{\bar{X}}_{n}=\mathbf{X}_{n}+\mathbf{\epsilon }</math>, for the next advection step,  achieving the final solution, <math>\mathbf{X}_{n+1}=A\left(\mathbf{\bar{X}}_{n}\right)</math>.    The steps are the following:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{l}  \mathbf{\hat{X}}_{n+1}=A\left(\mathbf{X}_{n}\right)=\mathbf{X}_{n}+\int_{t_{n}}^{t_{n+1}}\mathbf{V}\left(\mathbf{X}_{n},t\right)\mathrm{ d}t \\   \mathbf{\hat{X}}_{n}=A^{R}\left(\mathbf{\hat{X}}_{n+1}\right)=\mathbf{\hat{ X}}_{n+1}-\int_{t_{n}}^{t_{n+1}}\mathbf{V}\left(\mathbf{\hat{X}}_{n+1},t\right)\mathrm{d}t \\   \mathbf{\bar{X}}_{n}=\dfrac{\left(3\mathbf{X}_{n}-\mathbf{\hat{X}}_{n}\right)}{2} \\   \mathbf{X}_{n+1}=A\left(\mathbf{\bar{X}}_{n}\right)=\mathbf{\bar{X}}_{n}+\int_{t_{n}}^{t_{n+1}}\mathbf{V}\left(\mathbf{\bar{X}}_{n},t\right)  \mathrm{d}t \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
|}

where <math>A\left(\mathbf{\circ }\right)</math> involves the computation of the  velocity field of the tracer at its current position, interpolating from the  nodal values (Figure [[#img-32|32]]).

===7.3 RK4 method===

<div id='img-33'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-runge-kutta.png|600px|Fourth-order Runge-Kutta 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]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 33:''' Fourth-order Runge-Kutta 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
|}
<div id='img-34'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-zalesask.png|600px|Zalesask's benchmark after 10  revolutions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 34:''' Zalesask's benchmark after 10  revolutions.
|}
<div id='img-35'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-vel-quality.png|600px|Velocity effect on the weld  quality.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 35:''' Velocity effect on the weld  quality.
|}
<div id='img-36'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-fsw-particle.png|600px|Material particles after the weld with a threaded pin.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 36:''' Material particles after the weld with a threaded pin.
|}
According to the forth-order accurate RK4 method, the particle position at  time-step <math>n+1</math> is computed from the advection of the initial position by  four weighted incremental displacement at intermediate time-steps (Figure  [[#img-33|33]]):

<span id="eq-95"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \mathbf{X}_{n+1}=\mathbf{X}_{n}+\frac{1}{6}\left[\Delta \mathbf{X}^{\left( 1\right)}+2\Delta \mathbf{X}^{\left(2\right)}+2\Delta \mathbf{X}^{\left( 3\right)}+\Delta \mathbf{X}^{\left(4\right)}\right]   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
|}

where the incremental displacements are computed as

<span id="eq-96"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{l}  \Delta \mathbf{X}^{\left(1\right)}=\mathbf{V}\left(\mathbf{X}_{n},t_{n}\right)\Delta t \\   \Delta \mathbf{X}^{\left(2\right)}=\mathbf{V}\left(\mathbf{X}_{n}+\dfrac{1 }{2}\Delta \mathbf{X}^{\left(1\right)},t_{n}+\dfrac{\Delta t}{2}\right) \Delta t \\   \Delta \mathbf{X}^{\left(3\right)}=\mathbf{V}\left(\mathbf{X}_{n}+\dfrac{1 }{2}\Delta \mathbf{X}^{\left(2\right)},t_{n}+\dfrac{\Delta t}{2}\right) \Delta t \\   \Delta \mathbf{X}^{\left(4\right)}=\mathbf{V}\left(\mathbf{X}_{n}+\Delta   \mathbf{X}^{\left(3\right)},t_{n}+\Delta t\right)\Delta t \end{array}    </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (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 [[#img-34|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 [[#img-35|35]] demonstrates the velocity effect on the final weld  quality. It demonstrates that the rotational and advancing velocities can  not be selected arbitrarily.    Figure [[#img-36|36]] shows the final state of the material particles  (a colorful set) after FSW process.

==8 Residual stresses==

In this section, the two-step strategy for computing the residual stresses  is based upon combining global and local level analysis. The strategy is  shown in figure [[#img-37|37]].<div id='img-37'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-local-global.png|600px|Coupling strategy at local and global level.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 37:''' Coupling strategy at local and global level.
|}
The first step consists of computing the heat generated by the FSW tool at  the thermo-mechanically affected zone (TMAZ). This is performed within the  local level analysis.

===8.1 Step 1: 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 pin-shoulder) 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 [[#img-25|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 thermo-physical 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 non-linear behavior of the constitutive models that characterize the  material.    The output of the local analysis, which is the heat power (<math>\dot{P}</math>) to be  used at global level, is calculated once the steady-state is achieved. It is  computed taking into account the heat generated from both friction (<math> q_{fric} </math>) and plastic dissipation (<math>\dot{D}</math>) integrated over the TMAZ  domain:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \dot{P}=\int\nolimits _{{\Omega }_{TMAZ}}\dot{Q}\mathrm{d}V=\int\nolimits _{{\Omega }_{TMAZ}}\dot{D}\mathrm{d}V+\int\nolimits _{\partial {\Omega }_{TMAZ}}q_{fric}\mathrm{d}S  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (97)
|}

where <math>\dot{Q}</math> is the power density introduced into the system.    Thus, performing the local level analysis provides the input data for the  global study.

===8.2 Step 2: The global level analysis===

A simulation carried out at the 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 TMAZ at each time-step of the analysis. The  effects induced by the welding process on the structural behavior are the  target of the study. These are the distortions, residual stresses or  weaknesses along the welding line, among others. At global level, the  analysis data are the heat source, gravity load, clamping system and cooling  device (clamping, release and cooling down). In the model, the small  differences in the temperature resulting from the characteristics of the  retreating or advancing sides of the welds are not considered (Figure [[#img-38|38]]).    <div id='img-38'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-glob-heat.png|600px|Schematic of a global model with  a moving heat source]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 38:''' Schematic of a global model with  a moving heat source
|}
For the performance of the global level analysis, the definition of the  input power is required. Since no filler material is needed for joining in  FSW, the power input does not form part of the industrial setting and is  generally unknown. Apart from experiments and calibration, the power can be  approximated mathematically or by use of inverse analysis. In this work, it  is proposed to compute it from the local level simulation.    The principal effect of the material processing at the stirring zone on the  structure is the thermal expansion. The clamping system prevents the  expansion of the material due to the thermal deformations which transform  into plastic strains and, therefore, producing residual stresses in both  longitudinal and transversal directions in the weld and its vicinity.    The only data transferred between the two analysis levels is the heat  source, as the deformations are negligible in the size of stirring zone  comparing with the structure dimension.    Keeping the same input power obtained from local level analysis, the global  level FSW simulation can be performed for desired component dimensions and  welding paths, independently of the domain and mesh used in the local level  analysis.    <div id='img-39'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-tmaz.png|600px|TMAZ schematic]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 39:''' TMAZ schematic
|}
<div id='img-40'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-heat-source.png|600px|Schematic diagram of heat source model]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 40:''' Schematic diagram of heat source model
|}
The numerical simulation at global level needs an ad-hoc procedure in order  to apply the power energy to the elements representing the  Thermo-Mechanically Affected Zone (TMAZ) (Figure [[#img-39|39]]). Therefore, a  search algorithm is required to identify those elements at each time step of  the simulation. The total power input (<math>\dot{P}</math>) delivered within the  time-step, <math>\Delta t</math>, is distributed among the elements of the TMAZ volume  during the current time-step. The upper surface of the volume is a rectangle  with dimension of advancing speed (<math>v</math>) <math>\times </math> time step (<math>\Delta t</math>) <math> \times </math> pin diameter (<math>d_{p}</math>). The diameter of the profile of this volume  (<math>d_{0}</math>) decreases parabolically with depth <math>z</math> in the work-piece thickness  (<math>d</math>) (Figure [[#img-40|40]]):

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \begin{array}{cc}  d_{0}\left(z\right)=d_{p}\sqrt{\frac{z+d}{z}} & -d\leq z\leq 0 \end{array}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (98)
|}

The mesh resolution of the TMAZ in the global analysis needs to be  sufficient for the definition of the volume where the power input is  inserted. According to the search algorithm, the volumes (<math>V_{TMAZ}</math>) of all  the elements representing the TMAZ are added up as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  V_{TMAZ}=\sum \limits _{e=1}^{n_{e\in TMAZ}}V^{\left(e\right)}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (99)
|}

so that the heat source distribution is defined as (Figure [[#img-38|38]])

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>  \dot{Q}=\frac{\dot{P}}{V_{TMAZ}}  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (100)
|}

This power re-distribution preserves the total amount of energy input  independently of the mesh used (<span id='citeF-39'></span>[[#cite-39|[39]]]).<div id='img-41'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-local-2.png|600px|Mesh resolution at local (left) and global (right)  levels.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 41:''' Mesh resolution at local (left) and global (right)  levels.
|}
<div id='img-42'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-global-2.png|600px|Temperature field at local  (left) and global (right) levels.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 42:''' Temperature field at local  (left) and global (right) levels.
|}
<div id='img-43'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-comparison.png|600px|Temperature evolution obtained from local and global  level analyses compared with experiment]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 43:''' Temperature evolution obtained from local and global  level analyses compared with experiment
|}
<div id='img-44'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:draft_Nogué_795022754-residu-exp.png|600px|Comparison between the stresses  obtained from global level analysis and experiment at a line orthogonal to  the weld line.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 44:''' Comparison between the stresses  obtained from global level analysis and experiment at a line orthogonal to  the weld line.
|}

===8.3 Experimental validation===

In this section, firstly, the comparison between the temperature fields  obtained at local and global levels are shown to validate the local-global  coupling strategy. Then, the residual stresses computed using the strategy  are illustrated.    The meshes at local and global levels are shown in figure [[#img-41|41]]. The  total dissipated power is obtained at the steady state and used as the heat  input energy for the global level analysis. It is applied to the control  volume representing the TMAZ at each time-step of the analysis. The  temperature fields at local and global levels at the end of the process is  shown in figure [[#img-42|42]]. The comparison between the temperature  fields obtained from both local and global level analyses for different  lines parallel to the weld line at the bottom surface or thermocouples  placed at identical locations with those obtained from experiments are  presented in figure [[#img-43|43]].    The global level analysis provides the residual stress field. The comparison  between the variations of the residual stress (normal stresses in x  direction) along the transverse direction (along the line perpendicular to  the weld line) and experimental results are plotted in figure [[#img-44|44]] at the same location. Longitudinal stresses after fixture  release and cooling are depicted in this figure.

==9 Conclusion==

This paper presents the challenges and achievements in the numerical  simulation of coupled thermo-mechanical FSW problems. The existing  methodologies are designed to simulate FSW process either at global or local  level.    In the global level approach, the objective is to study the effect of a  moving heat source, in a Lagrangian framework, on the behavior of the whole  structure to be welded (the input energy is known).    In the local level approach, the objective is to analyze the heat affected  zone and to study the heat source locally. This can be used as energy input  for the global analysis (the heat generation comes from plastic dissipation  and frictional effects).    Generally, more attention has been paid to local level simulation as the  complex coupled thermo-mechanical phenomenon are mostly concentrated in the  HAZ. In the local level approach, the effect of high rotational velocity is  studied. Contact, strain-rate localization, strongly non-linear and  isochoric material behaviors and convective domination make the local  analysis more complicated and important.    In this paper, several aspects of FSW modeling are highlighted:

* '''Coupled thermo-mechnical model:''' A robust and accurate  thermo-mechanical solver for the simulation of FSW processes has been  developed. The transient thermal model proposed in both Lagrangian and  Eulerian kinematic frameworks allows for local and global level analysis.  Suitable boundary conditions based on Newton's law are introduced to deal  with heat conduction at contact interface and heat convection at the  external boundaries.
* '''Apropos kinematic framework:''' To study FSW process, a suitable  and robust kinematic setting has been proposed. The apropos kinematic  setting, combines ALE, Eulerian and Lagrangian frameworks for different  parts of the domain avoiding remeshing and variable remapping procedures.  Moreover, it permits to account for the arbitrary and complex pin shapes  typically found in FSW process.
* '''Mixed stabilized formulation for the mechanical problem:''' In  the hypothesis of isochoric material behavior as for liquid-like phase or  more generally when the deformations are mainly deviatoric (J2 plasticity  model used for metals), the stabilized mixed <math display="inline">\mathbf{v}/p</math> formulation has  been introduced. Due to the good performance, this technology must be used  especially for industrial simulations with triangular/tetrahedral meshes for  the domain discretization . Pressure locking is overcome using VMS  stabilization methods. Two stabilization methods are compared: ASGS and  OSGS. The OSGS stabilization method shows higher accuracy even if the  necessary computational effort is superior.
* '''Mixed stabilized formulation for convective dominated thermal  problem:''' To deal with thermal analysis where the convective term dominates  (ALE/Eulerian framework), the convection instabilities has been overcome  using VMS stabilization methods. This leads to a very accurate treatment of  this numerical problem. The stabilization method based on the sub-grid  scales approach circumvents the Ladyzenskaja-Babuska-Brezzi (LBB)  condition, allowing an accurate and robust formulation for linear elements.  Classical GLS and SUPG stabilization methods can be recovered as a  particular case of the sub-grid scale stabilization framework.
* '''Constitutive modeling:''' The thermo-elasto-visco-plastic  constitutive model used for the simulation of FSW process at global level  has been developed. This model covers the full temperature range leading to  an accurate evolution of temperature, distortions and the residual stresses.    Thermo-rigid-visco-plastic constitutive models such as Norton-Hoff, Carreau  or Sheppard-Wright have been used for the simulations of FSW processes at  local level. These kind of models neglect both the elastic and the thermal  strains. They are intended to solve mechanical problem dominated by very  high deviatoric strains as well as high strain rates. It is important to  observe that also in this case the thermal coupling is essential to capture  the real material behavior. This leads to highly non-linear models strongly  coupled with the temperature field, which can be treated as non-Newtonian  flows, with a temperature dependent viscosity.
* '''Strain localization problem: '''The FSW process produces a  strong strain localization in the stirring zone close to the pin. A sliding  zone concentrates in a thin layer under a pure shear stress field.    The solution of such problem is not feasible using a standard velocity-based  element. However, the use of mixed <math display="inline">\mathbf{v}/p</math> element is able to capture  the strain localization satisfactorily. Also in this case the proposed  VMS-based stabilization technique is suitable to treat such complex problem.
* '''Material flow visualization:''' The paper also deals with the  problem of material tracking in the stirring zone. The simulation is carried  out using a rigid visco-plastic constitutive model and this reduces into a  non-Newtonian fluid-like analysis where it is important to visualize the  material flow. The model provides this insight by tracking the particles  position during the full motion of the FSW process in the stirring zone.    The trend observed is consistent with the results of experiments, showing  the motion of the markers inserted into the weld. The model correctly  predicts the material flow around the retreating side of the pin. Results of  the simulation correctly demonstrate that the material is moved away from  the direction of travel and material sticking could occur at the tool  interface. The resulting material flow is not symmetric with respect to the  joint line and the flow patterns on the advancing and retreating sides are  different.    The model is able to include the effect of complex tool geometries as  threaded pins. Some distinctive 3D features of the flow can be appreciated.  These included: the material flow within the thickness direction and a  centrifugal material flow pattern around the pin.    Among the integration methods (BES, BFECC and RK4) used for the problem of  interest, RK4 is found to be the most accurate when a larger time-step is  used.
* '''Residual stress computation''': The paper proposes a  local-global two level strategy for the prediction of residual stresses in  FSW processes. At the first step, the heat power is calculated at local  level from both viscous dissipation and friction. Then the heat power is  transferred into the global level simulation in order to obtain the residual  stresses.

===BIBLIOGRAPHY===

<div id="cite-1"></div>
'''[[#citeF-1|[1]]]'''  Agelet de Saracibar, C., Cervera, M. and  Chiumenti, M. (1999) On the formulation of coupled thermoplastic problems  with phase-change, Int. J. of Plasticity 15:1&#8211;34.    <div id="cite-2"></div>
'''[[#citeF-2|[2]]]'''  Agelet de Saracibar, C., Chiumenti, M., Valverde, Q. and  Cervera, M. (2006) On the orthogonal sub-grid scale pressure stabilization  of finite deformation J2 plasticity, Comp. Method in Appl. Mech. and Eng.  195:1224-1251.    <div id="cite-3"></div>
'''[[#citeF-3|[3]]]'''  Agelet de Saracibar, C., Chiumenti, M., Santiago, D.,  Dialami, N. and Lombera, G. (2010) On the numerical modeling of fsw  processes. In Proceedings of the International Symposium on Plasticity and  its Current Applications, Plasticity 2010, St. Kitts, St. Kitts and Nevis,  January 3 - 8.    <div id="cite-4"></div>
'''[[#citeF-4|[4]]]'''  Agelet de Saracibar, C., Chiumenti, M., Santiago, D.,  Cervera, M., Dialami, N. and Lombera, G. (2010) A computational model for  the numerical simulation of fsw processes. In Proceedings of the 10th  International Conference on Numerical Methods in Industrial Forming  Processes 1252: 81 &#8211; 88, Pohang, South Korea, June 13-17. AIP Conference  Proceedings. doi: 10.1063/1.3457640.    <div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  Agelet de Saracibar, C., Chiumenti, M., Santiago, D.,  Cervera, M., Dialami, N. and Lombera, G. (2011) Advances in the numerical  simulation of 3D FSW processes. In Proceedings of the International  Symposium on Plasticity and its Current Applications, Plasticity 2011,  Puerto Vallarta, Mexico, January 3 - 8.    <div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  Agelet de Saracibar, C., Chiumenti, M., Cervera, M.,  Dialami, N. and Seret, A. (2014) Computational modeling and sub-grid scale  stabilization of incompressibility and convection in the numerical  simulation of friction stir welding processes. Archives of Computational  Methods in Engineering, 21(1), 3-37.    <div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  Agelet de Saracibar, C. and Chiumenti, M. (1998) On the  numerical modeling of frictional wear phenomena, Computer Methods in Applied  Mechanics and Engineering 177:401-426.    <div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  Askari, A., Silling, S., London, B. and Mahoney, M.  (2001) In 'Friction Stir Welding and Processing', (Ed. KV. Jata et al.)  43-54, Warrendale, PA, TMS.    <div id="cite-9"></div>
'''[[#citeF-9|[9]]]'''  Askari, A., Silling, S., London, B. and Mahoney, M.  (2003) modeling and Analysis of Friction Stir Welding Processes, The 4th  International Symposium on Friction Stir Welding, GKSS Workshop, Park City,  Utah, USA.    <div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  Assidi, M. and Fourment, L. (2009) Accurate 3D friction  stir welding simulation tool based on friction model calibration,  International Journal of Material Forming 2(1): 327-330.    <div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  Assidi, M., Fourment, L., Guerdoux, S. and Nelson, T.  (2010) Friction model for friction stir welding process simulation:  calibrations from welding experiments, Int. J. Mach. Tools Manuf.  50(2):143&#8211;155.    <div id="cite-12"></div>
'''[[#citeF-12|[12]]]'''  Bastier A, Maitournam MH, Roger F, Dang Van K. (2008)  Modelling of the residual state of friction stir welded plates. J Mater  Process Technol;200:25&#8211;37    <div id="cite-13"></div>
'''[[#citeF-13|[13]]]'''  Bendzsak, G., North, T. and Smith, C. (2000) An  Experimentally Validated 3D Model for Friction Stir Welding, Proceedings of  the 2nd International Symposium on Friction Stir Welding, Gothenburg, Sweden.    <div id="cite-14"></div>
'''[[#citeF-14|[14]]]'''  Bendzsak, G., North, T. and Smith, C. (2000) Material  Properties Relevant to 3-D FSW Modeling, Proceedings of the 2nd  International Symposium on Friction Stir Welding, Gothenburg, Sweden.    <div id="cite-15"></div>
'''[[#citeF-15|[15]]]'''  Bendzsak, O. J. and North, T. H. (1998) in  'Mathematical modeling of Weld Phenomena 4', (ed. H. C. Cerjak), 429-443,  London, 10M Communications Ltd.    <div id="cite-16"></div>
'''[[#citeF-16|[16]]]'''  Boman R., Ponthot J-P. (2012), Efficient ALE mesh  management for 3D quasi-Eulerian problems, International Journal for  Numerical Methods in Engineering, 92, 857-890,  http://dx.doi.org/10.1002/nme.4361    <div id="cite-17"></div>
'''[[#citeF-17|[17]]]'''  Boman R., Ponthot J-P (2013), Enhanced ALE data transfer  strategy for explicit and implicit thermomechanical simulations of  high-speed processes, International Journal of Impact Engineering, 53,  62-73, http://dx.doi.org/10.1016/j.ijimpeng.2012.08.007    <div id="cite-18"></div>
'''[[#citeF-18|[18]]]'''  Brezzi, F., and Fortin, M. (1991) Mixed and Hybrid Finite  Element Methods, Spinger, New York.    <div id="cite-19"></div>
'''[[#citeF-19|[19]]]'''  Buffa, G., Hu, J., Shivpuri, R. and Fratini, L. (2006) A  continuum based fem model for friction stir welding-model development,  Materials Science and Engineering A 419:389&#8211;396.    <div id="cite-20"></div>
'''[[#citeF-20|[20]]]'''  Buffa, G., Hu, J., Shivpuri, R. and Fratini, L. (2006)  Design of the friction stir welding tool using the continuum based FEM  model, Materials Science and Engineering A 419:381&#8211;388.    <div id="cite-21"></div>
'''[[#citeF-21|[21]]]'''  Buffa, G., Ducato, A. And Fratini L. (2011) Numerical  procedure for residual stresses prediction in friction stir welding, Finite  Elements in Analysis and Design 47:470&#8211;476.    <div id="cite-22"></div>
'''[[#citeF-22|[22]]]'''  Bussetta, P., Dialami, N., Boman, R., Chiumenti, M.,  Agelet de Saracibar, C., Cervera, M. and Ponthot, J.-P. (2013) Comparison of  a fluid and a solid approach for the numerical simulation of Friction Stir  Welding with a non-cylindrical pin, Steel research international,85 (6),  968-979.    <div id="cite-23"></div>
'''[[#citeF-23|[23]]]'''  Bussetta, P., Dialami, N., Chiumenti, M., Agelet de  Saracibar, C., Cervera, M. and Ponthot, J.-P. (2015), 3D numerical models of  FSW processes with non-cylindrical pin, AMPT 2014, 17th International  Conference on Advances in Material & Processing Technology, Atlantis The  Palm, Dubai, United Arab Emirates, 16-20 November 2014    <div id="cite-24"></div>
'''[[#citeF-24|[24]]]'''  Bussetta, P., Dialami, N., Chiumenti, M., Agelet de  Saracibar, C., Cervera, M., Boman, R. and Ponthot, J.-P. (2015) 3D numerical  models of FSW processes with non-cylindrical pin, Advanced Modeling and  Simulation in Engineering Sciences 2-27,  http://dx.doi.org/10.1186/s40323-015-0048-2, in press.    <div id="cite-25"></div>
'''[[#citeF-25|[25]]]'''  Bussetta P., Boman R., Ponthot J-P., Efficient 3D  data transfer operators based on numerical integration, Submitted to  International Journal for Numerical Methods in Engineering    <div id="cite-26"></div>
'''[[#citeF-26|[26]]]'''  Carreau, P.J. (1972) Rheological equations from  molecular network theories, Trans. Soc. Rheol. 16 (1):99&#8211;127.    <div id="cite-27"></div>
'''[[#citeF-27|[27]]]'''  Cederqvist, L. (2001) Properties of friction stir  welded aluminum lap joints. Master of Science Thesis, Department of  Mechanical Engineering, University of South Carolina, Columbia, USA.    <div id="cite-28"></div>
'''[[#citeF-28|[28]]]'''  Cervera, M., Agelet de Saracibar, C. and  Chiumenti, M. (1999) Thermo-mechanical analysis of industrial solidification  processes, Int. J. Numer. Meth. Eng. 46:1575&#8211;1591.    <div id="cite-29"></div>
'''[29]'''  Cervera, M., Chiumenti, M., Valverde, Q. and  Agelet de Saracibar, C. (2003) Mixed linear/linear simplicial elements for  incompressible elasticity and plasticity, Comp. Method in Appl. Mech. and  Eng. 192:5249-5263.    <div id="cite-30"></div>
'''[[#citeF-30|[30]]]'''  Cervera, M., Agelet de Saracibar, C. and Chiumenti, M.  (2002) Comet &#8211; a coupled mechanical and thermal analysis code. data input  manual. version 5.0. Technical Report IT-308, CIMNE, Barcelona, Spain. URL  http://www.cimne.com/comet.    <div id="cite-31"></div>
'''[[#citeF-31|[31]]]'''  Chao, Y.J., Qi, X. (1998) Thermal and  thermo-mechanical modeling of friction stir welding of aluminum alloy  6061-T6, J. Mater. Process. Manuf. Sci. 7:215&#8211;233.    <div id="cite-32"></div>
'''[[#citeF-32|[32]]]'''  Chao, Y.J., Qi, X. (1999) Heat transfer and  thermo-mechanical modeling of friction stir joining of AA6061-T6 plates,  Proceedings of the First International Symposium on Friction Stir Welding,  Thousand Oaks, CA, USA.    <div id="cite-33"></div>
'''[[#citeF-33|[33]]]'''  Chao, Y.J., Qi, X. and Tang, W. (2003) Heat  transfer in friction stir welding: experimental and numerical studies, ASME  J. Manuf. Sci. Eng. 125:138&#8211;145.    <div id="cite-34"></div>
'''[[#citeF-34|[34]]]'''  Chen, C. and Kovacevic, R. (2004) Thermomechanical modeling  and force analysis of friction stir welding by the finite element method,  Proc. Inst Mech. Eng., Part C: J. Mech. Eng. Sci. 218:509&#8211;519.    <div id="cite-35"></div>
'''[[#citeF-35|[35]]]'''  Chenot, J.L. and Massoni, E. (2006) Finite element  modeling and control of new metal forming processes, Journal of Machine  Tools and Manufacture, 46:1194&#8211;1200.    <div id="cite-36"></div>
'''[36]'''  Chiumenti, M., Valverde, Q., Agelet de  Saracibar, C. and Cervera, M. (2002) A stabilized formulation for elasticity  using linear displacement and pressure interpolations, Comp. Meth. in Appl.  Mech. and Eng. 191:5253-5264.    <div id="cite-37"></div>
'''[[#citeF-37|[37]]]'''  Chiumenti, M., Valverde, Q., Agelet de  Saracibar, C. and Cervera, M. (2004) A Stabilized Formulation for  Incompressible Plasticity Using Linear Triangles and Tetrahedra,  International Journal of Plasticity, 20:1487-1504.    <div id="cite-38"></div>
'''[[#citeF-38|[38]]]'''  Chiumenti, M., Agelet de Saracibar, C. and Cervera,  M. (2008) On the Numerical Modeling of the Thermo-mechanical Contact for  Casting Analysis, J. Heat Transfer 130:061301-1.    <div id="cite-39"></div>
'''[[#citeF-39|[39]]]'''  Chiumenti, M., Cervera, M., Salmi, A., Agelet de  Saracibar, C., Dialami, N., Matsui, K. (2010) Finite element modeling of  multi-pass welding and shaped metal deposition processes. Comput Method Appl  Mech Eng, 199:2343-2359.    <div id="cite-40"></div>
'''[[#citeF-40|[40]]]'''  Chiumenti, M., Cervera, M., Agelet de Saracibar, C.  and Dialami, N. (2012) Numerical modeling of friction stir welding  processes. Computer Methods in Applied Mechanics and Engineering,  254:353&#8211;369.    <div id="cite-41"></div>
'''[[#citeF-41|[41]]]'''  Chiumenti, M., Cervera, M., Agelet de Saracibar, C.  and Dialami, N. (2013) A novel stress accurate FE technology for highly  non-linear analysis with incompressibility constraint. application to the  numerical simulation of the FSWprocess. In Proceedings of the International  Conference on Numerical Methods in Forming Processes, NUMIFORM, Shenyang,  China. AIP Conference Proceedings.    <div id="cite-42"></div>
'''[42]'''  Cho, Y.I. and Kensey, K.R. (1991) Effects of the  non-Newtonian viscosity of blood on flows in a diseased arterial vessel.  Part I: Steady flows, Biorheology 28:241&#8211;262.    <div id="cite-43"></div>
'''[[#citeF-43|[43]]]'''  Cho, J., Boyce, D. and Dawson, P. (2005) modeling strain  hardening and texture evolution in friction stir welding of stainless steel,  Materials Science and Engineering A, 398:146&#8211;163.    <div id="cite-44"></div>
'''[[#citeF-44|[44]]]'''  Colegrove, P., Pinter, M., Graham, D. and  Miller, T. (2000) Three dimensional flow and thermal modeling of the  friction stir welding process, Proceedings of the Second International  Symposium on Friction Stir Welding, Gothenburg, Sweden, June 26&#8211;28.    <div id="cite-45"></div>
'''[[#citeF-45|[45]]]'''  Colegrove, P.A. (2004) modeling the friction stir  welding of aerospace alloys, 5th International Friction Stir Welding  Symposium, Metz, France, 14-16 September.    <div id="cite-46"></div>
'''[46]'''  Cross, M.M. (1965) Rheology of non-Newtonian fluids: a  new flow equation for pseudoplastic systems, J. Colloid Sci. 20:417&#8211;437.    <div id="cite-47"></div>
'''[[#citeF-47|[47]]]'''  Dang Van, K., Maitournam, M.H., (1993) Steady-state flow in  classical elastoplasticity: applications to repeated rolling and sliding  contact. J. Mech. Phys. Solids, 1691&#8211;1710    <div id="cite-48"></div>
'''[[#citeF-48|[48]]]'''  De Vuyst, T., D'Alvise, L., Robineau, A. and  Goussain, J.C. (2006) Material flow around a friction stir welding tool &#8211;  Experiment and Simulation, 8th International Seminar on Numerical Analysis  of Weldability, Graz, Austria.    <div id="cite-49"></div>
'''[[#citeF-49|[49]]]'''  De Vuyst, T., D'Alvise, L., Robineau, A. and  Goussain, J.C. (2006) Simulation of the material flow around a friction stir  welding tool, TWI 6th International Symposium on FSW, Saint-Sauveur.    <div id="cite-50"></div>
'''[[#citeF-50|[50]]]'''  Dialami, N., Chiumenti, M., Cervera, M. and Agelet de  Saracibar, C. (2013) An apropos kinematic framework for the numerical  modeling of friction stir welding. Computers and Structures 117:48&#8211;57.    <div id="cite-51"></div>
'''[[#citeF-51|[51]]]'''  Dialami, N., Chiumenti, M., Cervera, M., Agelet de  Saracibar, C. and Ponthot, J.-P. (2013) Material Flow Visualization in  Friction Stir Welding via Particle Tracing, International Journal of Metal  Forming 8 (2), 167-181.    <div id="cite-52"></div>
'''[[#citeF-52|[52]]]'''  Dialami, N., Chiumenti, M., Cervera, M., Agelet de  Saracibar, C., Ponthot, J.-P. and Bussetta, P. (2014) Numerical simulation  and visualization of material flow in Friction Stir Welding via Particle  Tracing, Numerical Simulations of Coupled Problems in Engineering, Springer  International Publishing,157-169.    <div id="cite-53"></div>
'''[[#citeF-53|[53]]]'''  Donea J., Huerta A., Ponthot J-P., Rodríguez-Ferr án A. (2004), Arbitrary Lagrangian-Eulerian Methods, in Encyclopedia of  Computational Mechanics, John Wiley & Sons, Ltd.,  http://dx.doi.org/10.1002/0470091355.ecm009    <div id="cite-54"></div>
'''[[#citeF-54|[54]]]'''  Dong, P., Lu, F., Hong, J. and Cao, Z. (2000) Analysis of  Weld Formation Process in Friction Stir Welding, The 2nd International  Symposium on Friction Stir Welding. Gothenburg, Sweden.    <div id="cite-55"></div>
'''[[#citeF-55|[55]]]'''  Dong, P., Lu, F., Hong, J.K., Cao, Z. (2001)  Coupled thermomechanical analysis of friction stir welding process using  simplified models, Sci. Technol. Weld. Jointing 6:281&#8211;287.    <div id="cite-56"></div>
'''[56]'''  Donne, C.D., Lima, E., Wegener,J., Pyzalla, A.  and Buslaps, T. (2001) Investigations on residual stresses in friction stir  welds, in: Proceedings of the Third International Symposium on Friction Stir  Welding, Kobe, Japan, September 27&#8211;28.    <div id="cite-57"></div>
'''[[#citeF-57|[57]]]'''  Ferrer, L.V. H. , Hernández, C. C. A. , Vargas  Aguilar, R. O. (2015) An Analytical Solution for Friction Stir Welding of an  AISI 1018 Steel. In Selected Topics of Computational and Experimental Fluid  Mechanics, Springer International Publishing, 473-480.    <div id="cite-58"></div>
'''[[#citeF-58|[58]]]'''  Fonda, R. W. and Lambrakos, S. G. (2002) Sci. Technol.  Weld. Joining. 7(3), 177-181.    <div id="cite-59"></div>
'''[[#citeF-59|[59]]]'''  Fourment, L., Guerdoux, S., Miles, M., Nelson, T.  (2004) Numerical Simulation of friction stir welding process using both  Lagrangian and Arbitrary Lagrangian Eulerian Formulations, 5 th  International Friction Stir Welding Symposium, Metz, France, 14-16 September.    <div id="cite-60"></div>
'''[[#citeF-60|[60]]]'''  Fourment, L. and Guerdoux, S. (2008) 3D numerical  simulation of the three stages of Friction Stir Welding based on friction  parameters calibration, International Journal of Material Forming 1(1):  1287-1290.    <div id="cite-61"></div>
'''[[#citeF-61|[61]]]'''  Fratini, L and Pasta, S. (2012) Residual stresses in  friction stir welded parts of complex geometry. Int J Adv Manuf Technol,  59:547&#8211;557.    <div id="cite-62"></div>
'''[[#citeF-62|[62]]]'''  Frigaard, O., Grong, O. and Midling, O.T.  (2001) A process model for friction stir welding of age hardening aluminum  alloys, Metal. Mater. Transition A 32:1189&#8211;1200.    <div id="cite-63"></div>
'''[[#citeF-63|[63]]]'''  Fu, L., Duan, L. Y. and Du, S. O. (2003) Numerical simulation  of inertia friction welding process by finite element method, Welding Joumal  82(3), 65s - 70s.    <div id="cite-64"></div>
'''[[#citeF-64|[64]]]'''  Ghanimi, Y., Celjak, H. C. and Faes, K. (2002) Modeling  of Friction Welding of Long Components, Proc. 6th Int. Conf. on 'Trends in  Welding Research', Phoenix, Arizona, ASM International, 329-333.    <div id="cite-65"></div>
'''[[#citeF-65|[65]]]'''  Gould, J. E. and Feng, Z. (1998) Heat flow model for  friction stir welding of aluminum alloys, Joumal of Materials Processing &  Manufacturing Science 7(2):185-194.    <div id="cite-66"></div>
'''[[#citeF-66|[66]]]'''  Grong, O. (1997) Metallurgical modeling of Welding, 2nd  edn, London, UK, Institute of Materials.    <div id="cite-67"></div>
'''[[#citeF-67|[67]]]'''  Guerdoux, S. and Fourment, L. (2009) A 3D numerical  simulation of different phases of friction stir welding, Modelling Simul.  Mater. Sci. Eng. 17 075001.    <div id="cite-68"></div>
'''[[#citeF-68|[68]]]'''  Hamilton, C., Dymek, S., Sommers, A. (2008) A thermal  model of friction stir welding in aluminium alloys, International Journal of  Machine Tools & Manufacture 48:1120-1130.    <div id="cite-69"></div>
'''[[#citeF-69|[69]]]'''  Heurtier, P., Jones, M.J., Desrayaud, C., Driver, J.H.,  Montheillet, F. and Allehaux, D. (2006) Mechanical and thermal modeling of  friction stir welding. J Mater Process Technol 171:348&#8211;357.    <div id="cite-70"></div>
'''[[#citeF-70|[70]]]'''  Hoff, H. J. (1954) Approximate Analysis of structures in  presence of moderately large creep deformation Quart. Appl. Math, 12:49.    <div id="cite-71"></div>
'''[[#citeF-71|[71]]]'''  Hughes, T.J.R. and Brooks, A.N. (1982) A theoretical  framework for Petrov-Galerkin methods, with discontinuous weighting  functions: applications to the streamline upwind procedure, in: R.H.  Gallagher, D.M. Norrie, J.T. Oden and O.C. Zienkiewicz, eds., Finite Element  in Fluids, vol. IV,Wiley, London 46-65..    <div id="cite-72"></div>
'''[[#citeF-72|[72]]]'''  Hughes, TJR., Franca, LP. and Hulbert, GM. (1989) A new  finite element formulation for computational fluid dynamics: VIII. The  Galerkin=least-squares method for advective&#8211;diffusive equations, Computer  Methods in Applied Mechanics and Engineering 73:173&#8211;189.    <div id="cite-73"></div>
'''[[#citeF-73|[73]]]'''  Hughes, T.J.R. (1995) Multiscale phenomena: Green's  function, Dirichlet-to Neumann formulation, sub-grid scale models, bubbles  and the origins of stabilized formulations, Comp. Meth. in Appl. Mech. and  Eng. 127:387-401.    <div id="cite-74"></div>
'''[[#citeF-74|[74]]]'''  Jorge Jr., A.M. and Balancín, O. (2004) Prediction of  Steel Flow Stresses under Hot Working Conditions, Materials Research  8(3):309-315.    <div id="cite-75"></div>
'''[[#citeF-75|[75]]]'''  Khandkar, M. Z. H. and Khan, J. A. (2001) Materials  Processing & Manufacturing Sci. 10(2):91-105.    <div id="cite-76"></div>
'''[[#citeF-76|[76]]]'''  Khandkar, M. Z. H., Khan, J. A. and Reynolds, A. P.  (2002) Proc. 6th Int. Conf. on 'Trends in Welding Research', Phoenix,  Arizona, ASM International, 218-223.    <div id="cite-77"></div>
'''[[#citeF-77|[77]]]'''  Khandkar, M. Z. H., Khan, J. A. and Reynolds, A. P.  (2003) Sci. Technol. Weld. Joining 8(3):166-174.    <div id="cite-78"></div>
'''[[#citeF-78|[78]]]'''  Khandkar, M. Z. H., Khan, J. A., Reynolds, A. P. and  Sutton, M.A. (2006) Predicting residual thermal stresses in friction stir  welded metals, Journal of Materials Processing Technology 174:195-203.    <div id="cite-79"></div>
'''[[#citeF-79|[79]]]'''  Lindner, K., Khandkar, Z., Khan, J., Tang, W. and  Reynolds, A. P. (2003) Proc. 4th Int. Symp. on 'Friction Stir Welding', Park  City, Utah, TWI Ltd.    <div id="cite-80"></div>
'''[[#citeF-80|[80]]]'''  Maol, A. and Massoni, E. (1995) Engineering Computations  12,497-512.    <div id="cite-81"></div>
'''[[#citeF-81|[81]]]'''  McClure, J. C., Feng, Z., Tang, T., Gould, J. E. and  Mun, L. E. (1998) Proc. 5th Int. Conf. on 'Trends in Welding Research', Pine  Mountain, Georgia, ASM International, 590-594.    <div id="cite-82"></div>
'''[[#citeF-82|[82]]]'''  McCune, R.W., Murphy, A., Price, M. and Butterfield, J.  (2012) The influence of friction stir welding process idealization on  residual stress and distortion predictions for future airframe assembly  simulations. J. Manuf. Sci. Eng. 134(3): 031011-031011-9.    <div id="cite-83"></div>
'''[[#citeF-83|[83]]]'''  Midling, O. T., and Grong, O. (1994) A process model for  friction welding of A1-Mg-Si alloys and Al-SiC metal matrix composites - I.  HAZ temperature and strain rate distribution, Acta Metall. Mater. 42(5):  1595&#8211;1609.    <div id="cite-84"></div>
'''[[#citeF-84|[84]]]'''  Midling, O.T. (1999) Effect of tool shoulder material  on heat input during friction stir welding. Proceedings of the First  International Symposium on Friction Stir Welding, Thousands Oaks, CA, USA,  June14&#8211;16.    <div id="cite-85"></div>
'''[[#citeF-85|[85]]]'''  Mitelea, I. and Radu, B. (1998) in 'Mathematical  modeling of Weld Phenomena 4', (ed. H. C. Cerjak), London, 10M  Communications Ltd, 444-454.    <div id="cite-86"></div>
'''[[#citeF-86|[86]]]'''  Nandan, R., Roy, G. G., Lienert, T. J. and DebRoy, T.  (2006) Numerical modeling of 3D plastic flow and heat transfer during  friction stir welding of stainless steel, Science and Technology of Welding  and Joining 11:526-537.    <div id="cite-87"></div>
'''[[#citeF-87|[87]]]'''  Nandan, R., Roy, G. G., Lienert, T. J. and DebRoy, T.  (2007) Three-dimensional heat and material flow during friction stir welding  of mild steel, Acta Materialia 55:883&#8211;895.    <div id="cite-88"></div>
'''[[#citeF-88|[88]]]'''  Nguyen, T. C. and Weckman, D. C. (1998) in 'Mathematical  modeling of Weld Phenomena 4', (ed. H. C. Celjak), London, 10M  Communications Ltd., 455-473.    <div id="cite-89"></div>
'''[[#citeF-89|[89]]]'''  Nikiforakis, N. (2005) Towards a whole system  simulation of FSW, 2nd FSW modeling and Flow Visualisation Seminar, GKSS  Forschungszentrum, Geesthacht.    <div id="cite-90"></div>
'''[[#citeF-90|[90]]]'''  Norton F.H. (1929) The creep of steel at hight  temperature. Mc Graw Hill, New York, USA.    <div id="cite-91"></div>
'''[[#citeF-91|[91]]]'''  Packer, S. M., Nelson, T. W. and Sorensen, C. D. (2001)  Tool and equipment requirements for friction stir welding ferrous and other  high melting temperature alloys. The 3rd International Symposium on Friction  Stir Welding, Kobe, Japan.    <div id="cite-92"></div>
'''[[#citeF-92|[92]]]'''  Paulo, RMF, Carlone, P, Valente, RAF, Teixeira-Dias, F  and Palazzo GS. (2013) Integrated design and numerical simulation of  stiffened panels including friction stir welding effects. Key Eng. Mater.  554&#8211;557:2237&#8211;42.    <div id="cite-93"></div>
'''[[#citeF-93|[93]]]'''  Pereyra, S., Lombera, G. A. and Urquiza, S. A. (2014).  Modelado numérico del proceso de soldadura FSW incorporando una té cnica de estimación de parámetros. Revista Internacional de Mé todos Numéricos para Cálculo y Diseño en Ingeniería, 30(3),  173-177.    <div id="cite-94"></div>
'''[[#citeF-94|[94]]]'''  Ponthot J.P. (2002), Unified stress update algorithms  for the numerical simulation of large deformation elasto-plastic and  elasto-viscoplastic processes, International Journal of Plasticity, 18,  91-126, http://dx.doi.org/10.1016/s0749-6419(00)00097-8    <div id="cite-95"></div>
'''[95]'''  Powell, R.E. and Eyring, H. (1944) Mechanism for  Relaxation Theory of Viscosity, Nature 154:427&#8211;428.    <div id="cite-96"></div>
'''[[#citeF-96|[96]]]'''  Rahmati Darvazi, A. andIranmanesh, M. (2014)  Prediction of asymmetric transient temperature and longitudinal . residual  stress in friction stir welding of 304L stainless steel, Materials &  Design, 55:812-820.    <div id="cite-97"></div>
'''[[#citeF-97|[97]]]'''  Reynolds, A.P., Khandkar, Z., Long, T. and Khan, J.  (2003) Utility of relatively simple models for understanding process  parameter effects on FSW, Source: Materials Science Forum  426-432(4):2959-2964.    <div id="cite-98"></div>
'''[[#citeF-98|[98]]]'''  Rosenthal, D. (1941) Mathematical theory of heat  distribution during welding and cutting, Welding J. 20 (5):220&#8211;234.    <div id="cite-99"></div>
'''[[#citeF-99|[99]]]'''  Russell, M. J. and Shercliff, H. R. (1998) Analytical  modeling of friction stir welding, Proc. 7th Int. Conf. on 'Joints in  Aluminium' - INALCO '98', Cambridge, UK, Abington Publishing, 197-207.    <div id="cite-100"></div>
'''[[#citeF-100|[100]]]'''  Russell, M.J. and Sheercliff, H.R. (1999)  Analytic modeling of microstructure development in friction stir welding.  Proceedings of the First International Symposium on Friction Stir Welding,  Thousands Oaks, CA, USA, June 14&#8211;16.    <div id="cite-101"></div>
'''[[#citeF-101|[101]]]'''  Santiago, D., Lombera, G., Urquiza, S., Agelet de  Saracibar, C. and Chiumenti, M. (2010) Modelado termo-mecánico del  proceso de friction stir welding utilizando la geometría real de la  herramienta. Revista Internacional de Métodos Numéricos para Cá lculo y Diseño en Ingeniería, 26:293 &#8211; 303.    <div id="cite-102"></div>
'''[[#citeF-102|[102]]]'''  Schmicker D., Naumenko K. and Strackeljan J. (2013) A  robust simulation of Direct Drive Friction Welding with a modified Carreau  fluid constitutive model, Computer Methods in Applied Mechanics and  Engineering , 265: 186-194.    <div id="cite-103"></div>
'''[[#citeF-103|[103]]]'''  Schmidt, H., Hattel, J. (2004) modeling thermo  mechanical conditions at the tool/matrix interface in Friction Stir welding,  The 5th International Symposium on Friction Stir Welding, Metz, France.    <div id="cite-104"></div>
'''[[#citeF-104|[104]]]'''  Seidel, T.U. and Reynolds, A.P. (2003) A 2-D  Friction Stir Welding Process Model Based on Fluid Mechanics, Sci. Technol.  Weld. Join. 8(3):175&#8211;183.    <div id="cite-105"></div>
'''[[#citeF-105|[105]]]'''  Sellars, CM. and Tegart, WJM. (1972) Hot workability,  Int. Metall. Rev. 17:1-24.    <div id="cite-106"></div>
'''[[#citeF-106|[106]]]'''  Sheppard, T. and Wright, D. (1979) Determination of  flow stress: Part 1 constitutive equation for aluminum alloys at elevated  temperatures, Met. Technol. 6:215-223.    <div id="cite-107"></div>
'''[[#citeF-107|[107]]]'''  Shi, Q-Y, Dickerson, T. and Shercliff, H. R. (2003) Heat  flow into friction stir tool, Proc. 4th Int. Symp. on 'Friction Stir  Welding', Park City, Utah, TWI Ltd. 108.    <div id="cite-108"></div>
'''[[#citeF-108|[108]]]'''  Song, M. and Kovacevic, R. (2002) A new heat transfer  model for friction stir welding, Trans. of NAMRI/SME XXX, 565-572.    <div id="cite-109"></div>
'''[[#citeF-109|[109]]]'''  Song, M., Kovacevic, R., Ouyang, J. and Valant, M. (2002)A  detailed three-dimensional transient heat transfer model for friction stir  welding, Proc. 6th Int. Conf. on 'Trends in Welding Research', Phoenix,  Arizona, ASM International, 212-217.    <div id="cite-110"></div>
'''[[#citeF-110|[110]]]'''  Song, M. and Kovacevic, R. (2003) Numerical and  experimental study of the heat transfer process in friction stir welding,  Proc. Instn Mech. Engrs. Part B: Journal of Engineering Manufacture  217:73-85.    <div id="cite-111"></div>
'''[[#citeF-111|[111]]]'''  Song, M. and Kovacevic, R. (2004) Heat transfer modeling  for both workpiece and tool in the friction stir welding process: A coupled  model, Proceedings of the Institution of Mechanical Engineers, Part B:  Journal of Engineering Manufacture 218(1):17-33.    <div id="cite-112"></div>
'''[[#citeF-112|[112]]]'''  Stewart, M. B., Nunes, A. C., Adams, G. P. and Romine,  P. (1998) A combined experimental and analytical modeling approach to  understanding friction stir-welding, Developments in Theoretical and Applied  Mechanics, XIX, 472.    <div id="cite-113"></div>
'''[[#citeF-113|[113]]]'''  Ulysse, P. (2002) Three-Dimensional modeling of the  Friction Stir Welding Process, International Journal of Machine Tools and  Manufacture 42:1549-1557.    <div id="cite-114"></div>
'''[114]'''  Walburn, F.J. and Schneck, D.J.(1976) A constitutive  equation for whole human blood, Biorheology 13:201&#8211;210.    <div id="cite-115"></div>
'''[[#citeF-115|[115]]]'''  Xu, S., Deng, X., Reynolds, A. and Seidel, T. (2001) Finite  Element Simulation of material flow in friction stir welding, Science and  Technology of Welding and Joining 6(3):191-193.    <div id="cite-116"></div>
'''[[#citeF-116|[116]]]'''  Xu, S. and Deng, X. (2003) Two and Three-Dimensional  Finite Element Models for the Friction Stir Welding Process, 4th  International Symposium on Friction Stir Welding, Park City, UT.    <div id="cite-117"></div>
'''[[#citeF-117|[117]]]'''  Xu, S. and Deng, X. (2004) Two and three-Dimensional  Finite Element models for the Friction Stir Welding Process, University of  South Carolina, Department of Mechanical Engineering, Columbia, SC 29208,  USA.    <div id="cite-118"></div>
'''[[#citeF-118|[118]]]'''  Yan, DY, Wu, AP, Silvanus, J. and Shi, QY (2011)  Predicting residual distortion of aluminum alloy stiffened sheet after  friction stir welding by numerical simulation. Materials & Design,  32(4):2284&#8211;2291.    <div id="cite-119"></div>
'''[119]'''  Yasuda, K., Armstrong, R.C. and Cohen, R.E. (1981)  Shear flow properties of concentrated solutions of linear and star branched  polystyrenes, Rheol. Acta 20:163&#8211;178.    <div id="cite-120"></div>
'''[[#citeF-120|[120]]]'''  Yoshikawa, K. (2003) A joining criterion for lap  joining of dissimilar metal materials of aluminum and stainless steel by  friction stir. The 4th International Symposium on Friction Stir Welding,  Park City, Utah, USA.    <div id="cite-121"></div>
'''[[#citeF-121|[121]]]'''  Zener, C. and Hollomon, J. H. (1944) Effect of strain rate  upon plastic flow of steel. Journal of Applied Physics 15:22&#8211;32.    <div id="cite-122"></div>
'''[[#citeF-122|[122]]]'''  Zhu, X.K. and Chao, Y.J. (2004). Numerical  simulation of transient temperature and residual stresses in friction stir  welding of 304L stainless steel, Journal of Materials Processing Technology  146: 263-272.