A finite element for laminar flow of incompressible fluids with inertia effects and thermomechanical coupling

The Friction Stir Spot Welding (FSSW) process involves large deformations in the neighborhood of the tool. The simulation of this process has to account for a pasty phase in which the material is stirred, and a phase remaining solid. An Arbitrary Lagrangian Eulerian (ALE) approach combined with respectively fluid and solid behaviours in each of those phases may allow to simulate a lot of rotations of the tool into the material while following the boundaries of the sheets. This work focuses on a first stage of this study, the development of a mixed formulation temperature/velocity/pressure of a fluid finite element P1+/P1 in the unsteady case.


Introduction
The Friction Stir Spot Welding process (Sakano et al., 2001) is derived from the Friction Stir Welding process (Thomas et al., 1991) which consists in creating a point weld between two sheets by penetration of a rotating tool into the material. Heat generated by friction and the motion of material lead to the generation of two phases, one in the neighborhood of the tool which becomes mushy and the other which remains solid. The unsteady nature of this process implies a motion of the interface between the two phases through phase change. The material mix generated ensures then the joint of parts after cooling.
To numerically model this process, standard lagrangian and eulerian approaches are not suitable. First, a completely lagrangian approach is not convenient because the material in the neighborhood of the tool undergoes large deformations, so one would have to frequently remesh resulting in the cumulation of errors and a huge numerical cost. Secondly, a usual eulerian approach would not allow to follow the boundaries of the sheets during the process.
In order to combine the advantages of both description, we plan to use an Arbitrary Eulerian Lagrangian (ALE) approach. The arbitrary motion will be set in such a way that the mesh will follow material points in the solid phase but not in the mushy zone. This will enable us to simulate important rotations of the tool while following the boundaries of the sheets, including the formation of a bead around the tool and the interface between both phases during the process.
Furthermore, in order to take into account the mushy state of the stirred material, this material will be described as a viscous fluid. Thus the model will couple an ALE numerical approach with two types of behaviour laws to correctly describe the fluid and the solid, including history effects in the solid.
We describe here the formulation of the P1+/P1 (Arnold et al., 1984) fluid finite element with thermomechanical coupling. It has already been implemented (Feulvarch et al., 2005;Feulvarch et al., 2007) in the code SYSWELD (2007) in a thermomechanical version for steady states, and is here extended to the unsteady case with a temperature/velocity/pressure formulation. The development of the fluid finite element described in section 2 represents the first step of the establishment of the ALE approach.
In section 3, verification of implementation of this element is made on the Couette viscosimeter problem. Two analytical solutions of this problem are developped respectively for an unsteady mechanical state, and for a coupling (unsteady thermal state) -(steady mechanical state). Each of these analytical solutions is developped with appropriate initial and boundary conditions in order to allow for easy comparisons with the numerical results. A good correlation is observed between numerical and analytical solutions on these test cases.
2. Mixed and unsteady temperature/velocity/pressure formulation for the P1+/P1 fluid finite element 2.1. The P1+/P1 finite element The modeling of the pasty phase as a viscous incompressible fluid involves the hydrostatic pressure as an unknown. From the numerical viewpoint, a particular discretization must be implemented to avoid some nonphysical numerical phenomena (Arnold et al., 1984;Hughes, 1987). The P1+/P1 tetrahedron is a linear finite element which ensures the continuity of the temperature, velocity and pressure fields. The "+" exponent refers to the "bubble" node added at the center of the element. This node allows to enrich the velocity field and to satisfy motion equations and the incompressibility condition. The velocity field is then interpolated on the element as: in which p is the current extremal node number, b refers to the bubble node and λ is the vector of degrees of freedom of the bubble node, homogeneous to velocities.

Weak form
The Friction Stir Spot Welding process must account for a strong thermomechanical coupling. The problem involves solving the equation of energy, the equations of motion and the internal incompressibility constraint. We shall assume that ∂Ω, the boundary of the domain Ω, admits the decomposition where ∂Ω v and ∂Ω F are respectively the parts of the boundary on which the velocities and tractions are prescribed. The temperature and the heat flux are, meanwhile, respectively prescribed on the parts ∂Ω θ and ∂Ω q of the boundary.
In the case of a mixed formulation, the weak form of the problem is stated as follows: given the body forces f , the heat source r, the prescribed tractions F d on ∂Ω F , and the thermal flux φ d on ∂Ω q , and initial conditions where the notation (•, ⋆) denotes the L 2 product of quantities • and ⋆, and other terms of the weak form are defined as follows The summation convention applies to repeated indices, and a comma stands for a partial differentiation (i.e., v i,j = ∂v i /∂x j ). The quantities D ij , k ij denote respectively the cartesian components of the strain rate tensor and the thermal conductivity tensor and D is the generated dissipation defined as where σ is the Cauchy stress tensor, related to the strain rate tensor through the constitutive law of the fluid: and p denotes the hydrostatic pressure, µ the kinematic viscosity and 1 the identity tensor. The trial solution and weighting spaces are defined as Notice that a stronger regularity is required for the temperature and velocity fields than for the pressure field.

Discretization
The weak form (W ) is spatially discretized using the finite element method. The resolution uses a backward Euler finite difference algorithm, the time interval is divided in N time steps ∆t n , n = 1 to N , T = N n=1 ∆t n . The implicit nature of the algorithm employed and the nonlinearities of this problem involve an iterative resolution. In the case of the Newton-Raphson method, the linearization of the semi-dicrete equations made at the k th iteration of the calculation leads to solve the following matricial system: n+1 are respectively the generalized mass, tangent stiffness and residue matrices, calculated at the k th iteration of the calculation. The residue and consistent stiffness matrices are defined as: and f int (k) n+1 are respectively the external and internal forces. These forces consist of a contribution of the q th node (i.e. q = 1, 4) of the tetrahedron and a contribution of the bubble node. The vector of the degrees of freedom q of the system is defined as: and its increment δq (k) as δq (k) = q The system of equations [9] is symbolically of the form: where K is composed of contributions of the sub-matrices: Thermal and mechanical convective terms are embedded in the residue vector. A special treatment of the thermal convective term ensures the stability of the numerical solution (Feulvarch et al., 2005;Bergheau et al., 2008;Belytschko et al., 2000).

Introduction of an extra approximation
In order to decrease the computational cost during the resolution, the degrees of freedom related to the bubble node are usually eliminated into equations related to extremal nodes. However, we can note that in the expression [13] we should solve ordinary differential equations on λ to make this elimination. So we introduce here an approximation which allows us to skip this further calculation.
The acceleration at the bubble node is given by: The relation 4 p=1 N (p) = 1 and the fact that the bubble node is located at the center of the tetrahedron element allow us to simplify the expression [15] using the fact that N (p) (x b ) = 1 4 . Moreover, the normalized shape function attached to the bubble node is such that N b (x b ) = 1. The acceleration expressed at the bubble node can thus be rewritten as: A usual approximation consists in lumping the mass matrix with a row sum technique. This approximation is valid when the acceleration field varies little on the element. We introduce here an analogous approximation but weaker, consisting in assuming that the acceleration at the bubble node is roughly equal to the average of those at the extremal nodes of this element. This leads to: Combining expressions [16] and [17] lead to set to zero the rates of the degrees of freedom of the bubble node: Notice we can also show that this approximation is compatible with a rigid body motion, the system of semi-discrete equations still verify the second Newton's law.
We can mention that an other approach of the formulation of this element has been made in (Basset, 2006), but the approximation made on the acceleration field introduced here seems clearer and more natural.
Finally, after temporal discretization, the resolution of the system is made with a linear system of equations of the type: 3. Verification of the fluid P1+/P1 finite element in unsteady conditions: the Couette viscosimeter The verification of the formulation of this element is made on the Couette viscosimeter problem. This problem is presented on figure 2. Two unsteady analytical solutions are respectively developped for an unsteady mechanical state, and for a coupling (unsteady thermal state) -(steady mechanical state). Solutions are obtained with appropriate initial and boundary conditions promoting the simplest expression of the solution fields, in order to allow for easy comparisons with the numerical results, rather than for the sake of physical interest.

Unsteady mechanical resolution
The assumptions of the unsteady mechanical solution are the following: -the problem is a transient axisymmetric case with a radial solution, -gravity is neglected, -the fluid is considered as Newtonian and the dynamic viscosity is fixed, independently of temperature and velocity, -the flow is incompressible, -the problem is isothermal, -the outer cylinder is fixed, the motion of the inner cylinder is imposed.
The calculation of the velocity field results from the tangential motion equation. Separating the variables, we obtain a Bessel equation (Carslaw et al., 1959), the solution of which is expressed as a linear combination of Bessel's functions of the first and second kinds at the first order in our case. Driving the inner cylinder velocity and adopting a convenient nonzero initial condition, this leads to the unsteady velocity field: where a and b are respectively the radii of the inner and outer cylinders, ρ and µ the density and the viscosity of the fluid, Ω 10 the initial rotation speed of the inner cylinder, J 1 and Y 1 the Bessel's functions of the first and second kind at the first order, and λ is a parameter homogeneous to the inverse of a length arising from the separation of variables, for which a value of unity (in the units used) will be chosen later.

Resolution for a coupling (unsteady thermal state) -(steady mechanical state)
In the equation of energy, the dissipation appears as a source term. This term is calculated from a problem chosen in such a way that the mechanical resolution be independent of the thermal solution, and can be made a priori. In our case, this term is given with the analytical solution of the purely mechanical steady state case. This combination leads to a coupled fluid/thermal problem in which only the thermal problem is transient. The thermal and mechanical constitutive parameters are independent of the temperature for this study.
With an appropriate change of the variable, it is possible to reduce the problem to a partial differential equation without a source term on a new variableT , the resolution is then made in the same manner as before. Dirichlet and Neumann homogeneous conditions are respectively imposed at r = a and r = b on the variableT . The unsteady thermal solution is then expressed as: where k, ρ, C are respectively the thermal conductivity, the density and the thermal capacity of the fluid, Ω 1 is the rotation speed of the inner cylinder and λ 1 the first root of the equation: which results from the boundary conditions. The homogeneous conditions applied on the variableT thus lead to non-homogeneous conditions by inverse variable change.

Comparison of numerical and analytical solutions
The fluid domain is discretized with P1+/P1 finite elements, a mesh is presented on figure 3. Inner and outer radii are fixed as a = 0.1m and b = 1m. For convenience, the material parameters of the fluid are all taken as unity.  The values of the analytical initial conditions are prescribed on each node, and velocities of the inner cylinder nodes are imposed according to an exponential. Thermal degrees of freedom have been blocked to a zero value at all nodes of the mesh. Figure  4 presents the comparison between the analytical and the numerical solution for the unsteady mechanical problem in which values at some time steps are plotted. The velocity field is plotted on a radial line, points identified on the plot correspond to finite element solutions extracted from the nodes of the mesh on a radial line apparent on the figure 3. We can observe a good correlation on intermediate time steps between both methods.

(Unsteady thermal state) -(steady mechanical state)
Initial thermal and mechanical conditions from the analytical solutions are prescribed on each node. Subsequently, the inner cylinder motion is driven at constant rotation speed Ω 1 = 10 rad/s, and temperature and flux conditions [21] are prescribed on nodes of the inner and outer cylinders. The figure 5 presents the comparison between analytical and numerical solutions for the unsteady thermal case. The temperature field is plotted on a radial line of the viscosimeter. We can again observe a good correlation between both solutions on the temperature field during the time steps of the calculation.

Conclusion
During the Friction Stir Spot Welding process, two phases are generated, one which is pasty and the other one solid. The work presented focuses on the development of the numerical tool designed for the simulation of the pasty phase. The P1+/P1 fluid finite element developped in the framework of a mixed formulation temperature/velocity/pressure is here extended to the unsteady case. The nonlinear finite element formulation is derived from the weak form of this strongly coupled problem. The time discretization is based on an implicit Euler finite difference scheme leading to an iterative resolution. To reduce the computational cost, the bubble node is eliminated before performing the whole resolution, through the introduction of an extra approximation on the acceleration field. This assumption allows us to avoid solving ordinary differential equations, while preserving the consistency of the semi-discrete equations with respect to a rigid body motion.
Afterwards, verification of this element is made on the Couette viscosimeter problem, two analytical solutions are developped to test the thermal and mechanical behaviours of this element. The comparison of both analytical and numerical solutions show a satisfactory behaviour of this element on these unsteady problems.