# Benchmark tests based on the Couette viscometer

Thomas Heuzé, Jean-baptiste Leblond, Jean-michel Bergheau

► To cite this version:

Thomas Heuzé, Jean-baptiste Leblond, Jean-michel Bergheau. Benchmark tests based on the Couette viscometer: I: Laminar flow of incompressible fluids with inertia effects and thermomechanical coupling. Computers & Mathematics with Applications, 2014, 67 (10), pp.1925-1937. 10.1016/j.camwa.2014.03.013. hal-01667116

HAL Id: hal-01667116

https://hal.science/hal-01667116

Submitted on 11 Oct 2018

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Benchmark tests based on the Couette viscometer—I: Laminar flow of incompressible fluids with inertia effects and thermomechanical coupling

Thomas Heuzé<sup>a,*</sup>, Jean-Baptiste Leblond<sup>b</sup>, Jean-Michel Bergheau<sup>c</sup>

The Couette viscometer is a well-known problem of fluid mechanics, well-suited for the verification of numerical methods. The aim of this work is to extend the classical steady state mechanical solution obtained in fluid mechanics, both to strongly-coupled thermomechanical problems in the case of laminar and incompressible fluid flows, and to solid-type nonlinear behaviours. Extended solutions will allow for the verification of new formulations of a mixed P1+/P1 finite element developed both in fluid and solid mechanics, within a temperature/velocity/pressure formulation coupled with an implicit (backward) Euler algorithm in time. In the present Part I, we address the case of the laminar flow of incompressible fluids with inertia effects and thermomechanical coupling. The verification performed on the reference solutions developed clearly evidence the good behaviour of the fluid finite element. The extension to solid-type nonlinear behaviours for strongly-coupled thermomechanical problems will be the subject of Part II.

## 1. Introduction

The Couette viscometer is a well-known apparatus designed to measure the viscosity of a fluid that consists of two coaxial rotating cylinders whose relative motion leads to shear the fluid enclosed between both cylinders. Though traditionally used in rheology, this viscometer is also a well-known problem of fluid mechanics and can be used as a benchmark test to verify the implementation of numerical methods. Its analytical steady solution for a laminar, isothermal and incompressible flow is reported in the early paper of Couette [1] in the case of a fixed inner cylinder and a rotating outer one, but was already addressed in the pioneering work of Stokes [2].

We intend in this work to extend the elementary solution in a simple way for numerical verification purposes. More specifically, we extend the solution to account for inertia effects and strong thermomechanical coupling. Two original analytical solutions are thus developed, for an unsteady mechanical state, and for a coupling (unsteady thermal state)–(steady

<footnote><sup>a</sup> Research Institute in Civil and Mechanical Engineering (GeM, UMR 6183 CNRS), École Centrale de Nantes, 1 rue de la Noë, F-44321 Nantes, France</footnote>

<footnote><sup>b</sup> Institut Jean Le Rond d'Alembert, Université Pierre et Marie Curie-Paris 6, UMR CNRS 7190, 4 Place Jussieu Tour 55-65, 75252 Paris, Cedex 05, France</footnote>

<footnote><sup>c</sup> Université de Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Etienne, Cedex 2, France</footnote>

<footnote>* Corresponding author. Tel.: +33 240372503. E-mail address: thomas.heuze@ec-nantes.fr (T. Heuzé).</footnote>
mechanical state) respectively, considering a Newtonian fluid.<fnref n="1" /> Each of these analytical solutions is developed with appropriate initial and boundary conditions in order to allow for easy comparisons with numerical results. These solutions are the topic of Section 2.

These new solutions permit to assess the implementation of an extension of the classical P1+/P1 finite element, or MINI element, initially introduced by Arnold et al. [5]. This element uses a mixed method which enables us to deal with the internal constraint of incompressibility. It has already been used in fluid mechanics [6], but also in the framework of material forming [7–10]. Moreover, it has already been implemented in fluid mechanics for a steady laminar flow, including a strong thermomechanical coupling [11–13] in the finite element code SYSWELD® [14], for welding applications. We present in this paper an extension of Feulvarch et al. [11,12]'s previous works on this element to the unsteady case in fluid mechanics, thus including convection terms, with a temperature/velocity/pressure formulation, coupled with an implicit (backward) Euler algorithm in time. The extension to transient states is performed by introducing a specific approximation of the acceleration field, to allow for the local elimination of the “bubble” degrees of freedom. It is shown that this approximation is compatible with an exact representation of rigid body motions. The formulation of this element is presented in Section 3.

The assessment of the fluid P1+/P1 finite element is subsequently performed using the reference solutions developed, and presented in Section 4. The aim of these solutions is to test the dynamic behaviour and the thermomechanical coupling embedded in the formulation of the finite element. Comparisons between analytical and numerical solutions are made on the transient velocity field for the first solution (unsteady mechanical case) and on the transient temperature field for the second solution (coupling (unsteady thermal state)–(steady mechanical state)). These comparisons show a good agreement between analytical and numerical solutions and evidence the good behaviour of the finite element.

Part II will be devoted to the extension of the Couette viscometer problem to solid-type nonlinear behaviours, in the case of strongly-coupled thermomechanical problems, and to the assessment of a solid P1+/P1 finite element developed within a (large displacement)–(large strain) framework, with a temperature/velocity/pressure formulation.

Such formulations appear to be desirable especially for some welding applications in which the material is stirred in the neighbourhood of a rotating tool, like in Friction Stir Welding-type processes. In such processes, the heat generated by friction and the motion of the material generates a mix which joins the parts to be welded after cooling. Thus, a strong interaction between thermal and mechanical effects drives the welding stage. Actually, the new elements are developed within the context of the modelling of the Friction Stir Spot Welding process [15,16] recently implemented in the finite element code SYSWELD® [14].

## 2. Analytical solutions of the Couette viscometer problem

Solutions developed in Sections 2.1 and 2.2 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. It should be noted that more physical initial and boundary conditions could be introduced, but would lead to solutions expressed as infinite sums, which proves less convenient for numerical verification.

We consider the Couette viscometer problem described on Fig. 1. The inner and outer cylinder radii are denoted a and b respectively. We shall consider in the following the outer cylinder as fixed and the motion of the inner one as driven.

### 2.1. Unsteady mechanical resolution

The laminar and incompressible flow of a Newtonian fluid is governed by the Navier–Stokes equations, written in the domain $\Omega$ as:

$$ \rho \left( \frac{\partial \mathbf{v}}{\partial t} + \nabla \mathbf{v} \cdot \mathbf{v} \right) = -\nabla p + \mathbf{f} + \mu \Delta \mathbf{v} \qquad (1) \\ \nabla \cdot \mathbf{v} = 0 $$

supplemented with appropriate initial and boundary conditions. In these equations $\rho$ and $\mu$ denote the density and the dynamic viscosity of the fluid, $p$, $\mathbf{f}$ and $\mathbf{v}$ the hydrostatic pressure, the body force and the velocity vector respectively. The symbols $\nabla$ and $\Delta$ denote the gradient and laplacian operators.

With regard to the Couette viscometer, the assumptions made in the unsteady mechanical solution are as follows:

* the problem is transient and axisymmetric,

* gravity forces are neglected,

* the dynamic viscosity is fixed, independently of temperature and velocity,

* the problem is isothermal.

<footnote n="1">¹ Unsteady mechanical solutions for the Couette problem have already been derived by Bernardin [3] and Kamran et al. [4], considering even more complex non-Newtonian fluid behaviours, but the solution expounded in Section 2.1 is expressed in a different and simpler format.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Parameterization of the Couette viscometer problem.</figcaption></figure>

The velocity field is orthonormal and of the form:

$$v(r) = v_{\theta}(r)e_{\theta} \qquad (2)$$

where $r \in [a, b]$ is the current radius. Also, the motion equations (1)$_1$ read:

$$\left\{ \begin{aligned} \frac{\rho v_{\theta}^2}{r} &= \frac{\partial p}{\partial r} \\ \rho \frac{\partial v_{\theta}}{\partial t} &= \mu \frac{\partial}{\partial r} \left( \frac{1}{r} \frac{\partial}{\partial r} (r v_{\theta}) \right). \end{aligned} \right. \qquad (3)$$

Eq. (3)$_2$ combined with appropriate initial and boundary conditions gives the velocity field, the hydrostatic pressure then derives from Eq. (3)$_1$. Assume now a velocity field in “separated form”:

$$v_{\theta}(r, t) = f(r)g(t). \qquad (4)$$

Introducing hypothesis (4) into Eq. (3)$_2$, one gets the following ordinary differential equations:

$$\begin{aligned} \frac{g'(t)}{g(t)} &= -\alpha \\ r^2 f''(r) + rf'(r) + \left(\frac{\rho\alpha}{\mu}r^2 - 1\right) f(r) &= 0 \end{aligned} \qquad (5)$$

where $\alpha$ is an arbitrary constant homogeneous to the inverse of a time. The differential equation (5)$_1$ is solved classically and gives an exponential solution fully determined by some initial condition. Eq. (5)$_2$ is a Bessel differential equation [17], whose solution is a linear combination of Bessel's functions of the first and second kinds. More precisely, we define the following change of variable and function:

$$x = \lambda r, \qquad F(x) = f(r) \qquad (6)$$

where $\lambda$ is homogeneous to the inverse of a length so that $x$ is dimensionless. Eq. (5)$_2$ can then be rewritten as:

$$x^2 F''(x) + x F'(x) + \left[ \frac{\rho \alpha}{\lambda^2 \mu} x^2 - 1 \right] F(x) = 0 \qquad (7)$$

which can be identified to be a Bessel equation of the first order by imposing the following relation between the constants $\alpha$ and $\lambda$:

$$\frac{\rho\alpha}{\lambda^2\mu} = 1. \qquad (8)$$

The solution of this equation is of the form:

$$F(x) = f(r) = C_1 J_1(\lambda r) + C_2 Y_1(\lambda r) \qquad (9)$$
where C₁ and C₂ are two constants to be determined from appropriate boundary conditions, J₁(λr) and Y₁(λr) being Bessel's functions of the first and second kinds at the first order.

Initial and boundary conditions are determined so as to be compatible with the assumption of separated variables. The boundary conditions of the problem are taken in the following form:

$$ 
\begin{aligned}
v_{\theta}(r=a, t) &= \Omega_0 a \exp(-\alpha t) \\
v_{\theta}(r=b, t) &= 0
\end{aligned} 
\qquad (10) 
$$

which implies by compatibility the following initial condition:

$$ v_{\theta}(r=a, t=0) = \Omega_0 a. \qquad (11) $$

The boundary conditions lead to the following system of equations:

$$ 
\begin{cases}
C_1 J_1(\lambda a) + C_2 Y_1(\lambda a) = \Omega_0 a \\
C_1 J_1(\lambda b) + C_2 Y_1(\lambda b) = 0.
\end{cases} 
\qquad (12) 
$$

From there follows the transient velocity field:

$$ v_{\theta}(r, t) = \Omega_0 a \frac{J_1(\lambda r) Y_1(\lambda b) - J_1(\lambda b) Y_1(\lambda r)}{J_1(\lambda a) Y_1(\lambda b) - J_1(\lambda b) Y_1(\lambda a)} \exp(-\alpha t) \quad \text{with } \alpha = \frac{\mu \lambda^2}{\rho} \quad (\lambda \text{ arbitrary parameter}). \qquad (13) $$

## 2.2. Resolution for a coupling (unsteady thermal state)–(steady mechanical state)

We consider now the Couette viscometer with a thermomechanical coupling. The thermal part of the solution is governed by the equation of energy conservation:

$$ \rho C \left( \frac{\partial T}{\partial t} + \nabla T \cdot \mathbf{v} \right) = \nabla \cdot (k \nabla T) + S \quad \forall \mathbf{x} \in \Omega \qquad (14) $$

supplemented with appropriate initial and boundary conditions. The symbols ρ, C and k denote the density, the heat capacity per unit mass and the thermal conductivity of the fluid, and S denotes a source term discussed below. This equation is combined with the Navier–Stokes equations (1) for the coupling. In this study, we assume that thermal and mechanical material parameters are fixed, independently of temperature and velocity. The kinematic assumptions introduced in the previous section are retained.

Accounting for the assumptions made, Eq. (14) can be simplified into:

$$ \rho C \frac{\partial T}{\partial t} = \frac{k}{r} \frac{\partial}{\partial r} \left( r \frac{\partial T}{\partial r} \right) + S. \qquad (15) $$

As already stated, the mechanical dissipation appears in Eqs. (14) and (15) as a source term S. We consider the case where this term may be computed from the mechanical solution alone; it then serves as an input in the thermal problem. More precisely, we assume this term to be given by the solution of the steady state mechanical problem and thus be time-independent.

These assumptions lead to a coupled fluid/thermal problem in which the sole thermal problem is transient.

The dissipation generated reads:

$$ S = 2\mu\mathbf{D} : \mathbf{D} \qquad (16) $$

where μ is the dynamic viscosity of the fluid and D the eulerian strain rate tensor. The velocity field, in the steady state case, reads:

$$ v_{\theta}(r) = \frac{\Omega_1 a^2}{b^2 - a^2} \left( \frac{b^2}{r} - r \right) \qquad (17) $$

where Ω₁ denotes the angular velocity imposed on the inner cylinder. This allows us to express the mechanical dissipation as:

$$ S(r) = \frac{A}{r^4}, \quad \text{with} \quad A = 4\mu\Omega_1^2 \frac{(ab)^4}{(b^2 - a^2)^2}. \qquad (18) $$

The main difficulty encountered in the solution of Eq. (15) is the presence of the source term which prevents the use of a technique of separation of variables. However, it is possible to reduce the problem to a partial differential equation without any source term on a new variable T, defined by:

$$ \bar{T}(r, t) = T(r, t) + h(r) \qquad (19) $$
where h(r) is a function to be defined. Introducing the decomposition (19) into Eq. (15), it is quite easy to identify h(r) so as to get a partial differential equation without any source term; it suffices to take:

$$h(r) = \frac{A}{4kr^2}. \tag{20}$$

The thermal problem can then be rewritten on the variable $\bar{T}$ as:

$$\rho C \frac{\partial \bar{T}}{\partial t} = \frac{k}{r} \frac{\partial}{\partial r} \left( r \frac{\partial \bar{T}}{\partial r} \right). \tag{21}$$

The technique of separation of variables is now applicable, the resolution is straightforward. Set

$$\bar{T}(r, t) = f(r)g(t). \tag{22}$$

Introducing hypothesis (22) into Eq. (21) leads to the following ordinary differential equations:

$$\begin{gathered} \frac{g'(t)}{g(t)} = -\alpha \\ f''(r) + \frac{f'(r)}{r} + \frac{\alpha}{\kappa} f(r) = 0 \end{gathered} \tag{23}$$

where $\alpha$ is an arbitrary constant, and $\kappa$ the thermal diffusivity defined as:

$$\kappa = \frac{k}{\rho C}. \tag{24}$$

One may observe that Eq. (23)$_2$ can be identified to be a Bessel equation, but this time at the order zero. The identification is made using the change of variable and function defined by Eq. (6), which leads to:

$$x^2 F''(x) + x F'(x) + \frac{\alpha x^2}{\kappa \lambda^2} F(x) = 0. \tag{25}$$

This equation is a Bessel equation of order zero provided the constants $\alpha$ and $\lambda$ are related by:

$$\frac{\alpha}{\kappa \lambda^2} = 1. \tag{26}$$

The solution reads:

$$F(x) = C_1 J_0(x) + C_2 Y_0(x) \tag{27}$$

where $C_1$ and $C_2$ are two constants to be determined from appropriate boundary conditions, $J_0(x)$ and $Y_0(x)$ being Bessel's functions of the first and second kinds at order zero.

Boundary conditions are defined so as to lead to the simplest possible expression of the solution field. Homogeneous Dirichlet conditions on $\bar{T}$ are therefore prescribed on both cylinders of the viscometer:

$$\begin{aligned} \bar{T}(r=a, t=0) &\Rightarrow f(a)=0 &\Rightarrow F(\lambda a) &= 0 \\ \bar{T}(r=b, t=0) &\Rightarrow f(b)=0 &\Rightarrow F(\lambda b) &= 0. \end{aligned} \tag{28}$$

These boundary conditions give nonzero imposed values when going back to the temperature $T$:

$$\begin{aligned} T(a, t) &= -\frac{A}{4ka^2} \\ T(b, t) &= -\frac{A}{4kb^2}. \end{aligned} \tag{29}$$

Condition (28) leads to the following system of equations:

$$\begin{bmatrix} J_0(\lambda a) & Y_0(\lambda a) \\ J_0(\lambda b) & Y_0(\lambda b) \end{bmatrix} \begin{pmatrix} C_1 \\ C_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}. \tag{30}$$

The determinant of the matrix on the left-hand side must vanish for a non-trivial solution to exist:

$$z(\lambda) = J_0(\lambda a) Y_0(\lambda b) - J_0(\lambda b) Y_0(\lambda a) = 0. \tag{31}$$

The constants $C_1$ and $C_2$ depend on the roots $\lambda_i$ of Eq. (31), which cannot be derived explicitly. Therefore, the roots of (31) have to be computed numerically for given values of the inner and outer radii of the viscometer, the solution becomes therefore semi-analytical. In the sequel, the inner radius is set to $a = 0.1$ m and the outer radius is set to $b = 1$ m. With these numerical values, the function $z(\lambda)$ is plotted in Fig. 2.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Plot of the function z(λ) on the interval λb ∈ [0, 20].</figcaption></figure>

Combining Eqs. (30) and (31), one gets:

$$ \eta = \frac{C_1}{C_2} = -\frac{Y_0(\lambda_i a)}{J_0(\lambda_i a)} = -\frac{Y_0(\lambda_i b)}{J_0(\lambda_i b)} \quad (32) $$

where λᵢ denotes the chosen root of equations (30). Therefore, the spatial form of the solution reads:

$$ \bar{T}(r, t) = K \left( J_0(\lambda_i r) + \frac{Y_0(\lambda_i r)}{\eta} \right) \exp(-\alpha t) \quad \text{with } \alpha = \kappa \lambda_i^2 \quad (33) $$

where K is a constant homogeneous to a temperature, arising from the integration of Eq. (23)₁. Once the solution T̄ is known, we can come back to the original problem posed on the temperature T by means of the inverse change of variable defined by (19). The transient temperature field of the thermomechanically coupled viscometer problem is expressed as:

$$ \begin{gathered} T(r, t) = K \left( J_0(\lambda_i r) + \frac{Y_0(\lambda_i r)}{\eta} \right) \exp(-\alpha t) - \frac{A}{4kr^2} \\ \text{with } \eta = -\frac{Y_0(\lambda_i a)}{J_0(\lambda_i a)}, \quad A = 4\mu\Omega_1^2 \frac{(ab)^4}{(a^2 - b^2)^2} \alpha = \kappa\lambda_i^2, \quad \kappa = \frac{k}{\rho C} \end{gathered} \quad (34) \\ T(a, t) = -\frac{A}{4ka^2} \\ T(b, t) = -\frac{A}{4kb^2} \\ T(r, t=0) = T_0(r). \end{gathered} $$

Comments

* Homogeneous Dirichlet boundary conditions (Eqs. (28)) have been imposed on T̄, this has led to the homogeneous equations (30) and thus to an equation governing conditions of existence of a non-trivial solution (31). It is possible to solve this problem with other couples of homogeneous boundary conditions. If we seek for a simple expression of the solution field, we can prescribe a combination of zero temperature and flux on the cylinders or only zero fluxes on both cylinders. In these cases, it turns out that the form of the analytical solution (Eq. (34)) remains the same, the change in the solution occurs only through a change of Eq. (31) defining the possible values of λ.

* Non-homogeneous boundary conditions imposed on T̄ would lead to a non-homogeneous system of equations in place of (30) that would admit a unique solution for constants C₁ and C₂, for all values of λ leading to a nonzero determinant.

3. Mixed and unsteady temperature/velocity/pressure formulation for the P1+/P1 fluid finite element

3.1. Weak form of the equations

In welding applications, a strong thermomechanical coupling occurs due to the importance of the mechanical dissipation. Therefore, one must simultaneously solve the equations of conservation of energy and momentum, with the internal
constraint of incompressibility. We shall assume that the boundary $\partial\Omega$ of the domain $\Omega$ admits the decompositions

$$ \partial\Omega = \partial\Omega_v \cup \partial\Omega_F \qquad \emptyset = \partial\Omega_v \cap \partial\Omega_F \qquad (35) $$

$$ \partial\Omega = \partial\Omega_{\theta} \cup \partial\Omega_{q} \qquad \emptyset = \partial\Omega_{\theta} \cap \partial\Omega_{q} \qquad (36) $$

where $\partial\Omega_v$ and $\partial\Omega_F$ denote respectively these parts on which the velocities and tractions are prescribed, and $\partial\Omega_{\theta}$ and $\partial\Omega_q$ those on which the temperature and heat fluxes are prescribed.

In a mixed formulation, the weak form of the problem reads as follows: given the body forces $\mathbf{f}$, the heat source $r$, the prescribed tractions $\mathbf{F}_d$ on $\partial\Omega_F$, and the thermal flux $\phi^d$ on $\partial\Omega_q$, plus some initial conditions $T_0(\mathbf{x})$, $\mathbf{v}_0(\mathbf{x})$,

$$ (W) \quad \left\{ \begin{array}{l} \text{Find}(\mathbf{T}, \mathbf{v}, \mathbf{p}) \in (\mathcal{T}_{ad} \times \mathcal{V}_{ad} \times \mathcal{P}_{ad}), \quad \forall t \in [0, T], \text{ such that} \\ \quad \forall (T^*, \mathbf{v}^*, \mathbf{p}^*) \in (\mathcal{T}_{ad}^0 \times \mathcal{V}_{ad}^0 \times \mathcal{P}_{ad}), \\ \quad \begin{aligned} & -\int_{\Omega} k \nabla T \cdot \nabla T^* d\Omega + \int_{\partial\Omega_q} \phi^d T^* dS + \int_{\Omega} 2\mu \mathbf{D} : \mathbf{D} T^* d\Omega + \int_{\Omega} r T^* d\Omega \\ &= \int_{\Omega} \rho C \frac{\partial T}{\partial t} T^* d\Omega + \int_{\Omega} \rho C (\nabla T) \cdot \mathbf{v} T^* d\Omega \\ &-\int_{\Omega} 2\mu \mathbf{D} : \mathbf{D}^* d\Omega + \int_{\Omega} p \nabla \cdot \mathbf{v}^* d\Omega + \int_{\partial\Omega_F} \mathbf{F}_d \cdot \mathbf{v}^* dS + \int_{\Omega} \mathbf{f} \cdot \mathbf{v}^* d\Omega \end{aligned} \\ \quad = \int_{\Omega} \rho \frac{\partial \mathbf{v}}{\partial t} \cdot \mathbf{v}^* d\Omega + \int_{\Omega} \rho (\nabla \mathbf{v}) \cdot \mathbf{v} \cdot \mathbf{v}^* d\Omega \\ \quad \int_{\Omega} p^* \nabla \cdot \mathbf{v} d\Omega = 0 \\ \quad T(\mathbf{x}, t=0) = T_0(\mathbf{x}) \\ \quad \mathbf{v}(\mathbf{x}, t=0) = \mathbf{v}_0(\mathbf{x}) \end{array} \right. \qquad (37) $$

where $\mathcal{T}_{ad}$, $\mathcal{V}_{ad}$ and $\mathcal{P}_{ad}$ are temperature, velocity and pressure function spaces, while $\mathcal{T}_{ad}^0$ and $\mathcal{V}_{ad}^0$ refer to the weighting spaces associated, that is with homogeneous Dirichlet boundary conditions.

## 3.2. The P1+/P1 finite element

The P1+/P1 finite element initially proposed by Arnold et al. [5] ensures in the proposed version the continuity of temperature, velocity and pressure fields. The velocity field approximation is enriched with a bubble function that enables this element to satisfy the Brezzi-Babuska conditions [5], and is thus interpolated as:

$$ \mathbf{v}^h(\mathbf{x}) = N^{(p)}(\xi, \eta, \zeta)\mathbf{v}^{(p)} + N^{(b)}(\xi, \eta, \zeta)\lambda \quad 1 \le p \le 4 \qquad (38) $$

where $N^{(p)}$ and $N^{(b)}$ denote the shape functions associated respectively to the current vertex node $p$, and the bubble node $b$; $\lambda$ is the vector of degrees of freedom of the bubble node, homogeneous to some velocity. Implicit summation is used in Eq. (38) with respect to the index $p$.

## 3.3. Semidiscrete equations

The weak form (37) is discretized with the P1+/P1 finite element; this leads to the following system of semidiscrete equations:

$$ \mathbf{M}\dot{\mathbf{q}} + \mathbf{f}^{\text{conv}} + \mathbf{f}^{\text{int}} = \mathbf{f}^{\text{ext}} \qquad (39) $$

where $\mathbf{q}$ is the vector of degrees of freedom of the system defined as:

$$ \mathbf{q}^T = \{\mathbf{T} \ \mathbf{v} \ \mathbf{p} \ \lambda\}. \qquad (40) $$

The internal forces $\mathbf{f}^{\text{int}}$, external forces $\mathbf{f}^{\text{ext}}$, convection forces $\mathbf{f}^{\text{conv}}$ are assembled from elementary quantities:

$$ \mathbf{f}^{\text{int}} = \sum_{e=1}^{N_e} \left\{ \begin{Bmatrix} \int_{T_e}^{\text{int}, (p)} \\ \int_{v_e}^{\text{int}, (p)} \\ \int_{p_e}^{\text{int}, (p)} \end{Bmatrix} \right\}, \quad \mathbf{f}^{\text{ext}} = \sum_{e=1}^{N_e} \left\{ \begin{Bmatrix} \int_{T_e}^{\text{ext}, (p)} \\ \int_{v_e}^{\text{ext}, (p)} \\ 0 \\ \int_{p_e}^{\text{ext}, (p)} \end{Bmatrix} \right\}, \quad \mathbf{f}^{\text{conv}} = \sum_{e=1}^{N_e} \left\{ \begin{Bmatrix} \int_{T_e}^{\text{conv}, (p)} \\ \int_{v_e}^{\text{conv}, (p)} \\ 0 \\ \int_{p_e}^{\text{conv}, (p)} \end{Bmatrix} \right\} \qquad (41) $$
and consist of components pertaining to temperature, velocity and pressure equations, associated to the degrees of freedom of the vertex nodes (denoted p, 1 ≤ p ≤ 4), plus components relative to bubble velocity equations associated to degrees of freedom of the bubble node (denoted b). The components of internal forces are detailed below:

$$
\begin{align}
f_T^{\text{int},(p)} &= - \int_{\Omega_e} \nabla N^{(p)} \cdot \mathbf{h} \, d\Omega \nonumber \\
\mathbf{f}_v^{\text{int},(p)} &= \int_{\Omega_e} 2\mu \mathbf{D} : \mathbf{D}(N^{(p)}\mathbf{e}_i) d\Omega - \int_{\Omega_e} p \nabla N^{(p)} d\Omega \nonumber \\
&= \left| \begin{aligned}[t]
& \left[ \left( \int_{\Omega_e} \mu \nabla N^{(q)} \otimes \nabla N^{(p)} d\Omega \right) + \left( \int_{\Omega_e} \mu \nabla N^{(q)} \cdot \nabla N^{(p)} d\Omega \right) \mathbf{1} \right] \cdot \mathbf{v}^{(q)} \\
& + \left[ \left( \int_{\Omega_e} \mu \nabla N^{(b)} \otimes \nabla N^{(p)} d\Omega \right) + \left( \int_{\Omega_e} \mu \nabla N^{(b)} \cdot \nabla N^{(p)} d\Omega \right) \mathbf{1} \right] \cdot \lambda \\
& - \int_{\Omega_e} p \nabla N^{(p)} d\Omega
\end{aligned} \right| \tag{42} \\
f_p^{\text{int},(p)} &= - \int_{\Omega_e} N^{(p)} \nabla \cdot \mathbf{v} \, d\Omega \nonumber \\
\mathbf{f}_b^{\text{int}} &= \int_{\Omega_e} 2\mu \mathbf{D} : \mathbf{D}(N^{(b)}\mathbf{e}_i) d\Omega - \int_{\Omega_e} p \nabla N^{(b)} d\Omega \nonumber \\
&= \left| \begin{aligned}[t]
& \left[ \left( \int_{\Omega_e} \mu \nabla N^{(q)} \otimes \nabla N^{(b)} d\Omega \right) + \left( \int_{\Omega_e} \mu \nabla N^{(q)} \cdot \nabla N^{(b)} d\Omega \right) \mathbf{1} \right] \cdot \mathbf{v}^{(q)} \\
& + \left[ \left( \int_{\Omega_e} \mu \nabla N^{(b)} \otimes \nabla N^{(b)} d\Omega \right) + \left( \int_{\Omega_e} \mu \nabla N^{(b)} \cdot \nabla N^{(b)} d\Omega \right) \mathbf{1} \right] \cdot \lambda \\
& - \int_{\Omega_e} p \nabla N^{(b)} d\Omega
\end{aligned} \right| \nonumber
\end{align}
$$

where $\Omega_e$ is the volume of the element, **1** the second-rank identity tensor, and **h** the heat flux vector linked to the temperature through the Fourier's law **h** = $-k\nabla T$, $k$ being the thermal conductivity. External forces are given by:

$$
\begin{align}
f_T^{\text{ext},(p)} &= \int_{\partial \Omega_e \cap \partial \Omega_q} \phi^d N^{(p)} dS + \int_{\Omega_e} r \bar{N}^{(p)} d\Omega + \int_{\Omega_e} (2\mu \mathbf{D} : \mathbf{D}) \bar{N}^{(p)} d\Omega \nonumber \\
\mathbf{f}_v^{\text{ext},(p)} &= \int_{\partial \Omega_e \cap \partial \Omega_F} \mathbf{F}_d N^{(p)} dS + \int_{\Omega_e} \mathbf{f} N^{(p)} d\Omega \tag{43} \\
\mathbf{f}_b^{\text{ext}} &= \int_{\Omega_e} \mathbf{f} N^{(b)} d\Omega. \nonumber
\end{align}
$$

Notice that $\int_{\partial\Omega_e \cap \partial\Omega_F} \mathbf{F}_d N^{(b)} dS = 0$ since the bubble function vanishes on the element boundary. The function $\bar{N}^{(p)}$ introduced in (43) is associated to the Petrov-Galerkin method [18-20] used on the thermal equations to avoid numerical instabilities as the Peclet number increases. This function is chosen in the case of an incompressible flow as:

$$
\bar{N}^{(p)} = N^{(p)} + \tau \mathbf{v} \cdot \nabla N^{(p)} \quad (44)
$$

where $\tau$ is an adjustable parameter. The same type of numerical instabilities is expected for the mechanical equations, governed by the Reynolds number. However in welding applications, the flow is pasty and the Reynolds number is small. Therefore, the Petrov-Galerkin method is not used for the mechanical equations. Finally, convection forces read:

$$
\begin{align}
f_T^{\text{conv},(p)} &= \int_{\Omega_e} \rho C(\nabla T) \cdot \mathbf{v} \bar{N}^{(p)} d\Omega \nonumber \\
\mathbf{f}_v^{\text{conv},(p)} &= \int_{\Omega_e} \rho (\nabla \mathbf{v}) \cdot \mathbf{v} N^{(p)} d\Omega \nonumber \\
\mathbf{f}_b^{\text{conv}} &= \int_{\Omega_e} \rho (\nabla \mathbf{v}) \cdot \mathbf{v} N^{(b)} d\Omega. \tag{45}
\end{align}
$$

Convection and internal forces may be expressed using a convective matrix **L** and a stiffness matrix **K** as:

$$
\left\{
\begin{array}{l}
\mathbf{f}^{\text{conv}} = \mathbf{Lq} \\
\mathbf{f}^{\text{int}} = \mathbf{Kq}
\end{array}
\right.
\qquad (46)
$$
where **q** is the vector of degrees of freedom of the system already defined in (40). The system of semidiscrete equations (39) may be written more explicitly, gathering matrices **L** and **K**:

$$
\sum_{e=1}^{N_e} \left( \begin{bmatrix} M_T^{pq} & \mathbf{0} & 0 & \mathbf{0} \\ \mathbf{0} & M_v^{pq} & \mathbf{1} & 0 \\ 0 & \mathbf{0} & 0 & \mathbf{0} \\ \mathbf{0} & M_v^{bq} & \mathbf{1} & 0 \end{bmatrix} \begin{bmatrix} \dot{T}^{(q)} \\ \dot{\mathbf{v}}^{(q)} \\ \dot{p}^{(q)} \\ \dot{\lambda} \end{bmatrix} + \begin{bmatrix} L_T^{pq} + K_T^{pq} & \mathbf{0} & 0 & \mathbf{0} \\ \mathbf{0} & L_v^{pq} \mathbf{1} + \mathbf{K}_v^{pq} & \mathbf{k}_v^{pq} & L_v^{pb} \mathbf{1} + \mathbf{K}_v^{pb} \\ 0 & (\mathbf{k}_p^{pq})^T & 0 & (\mathbf{k}_p^{pb})^T \\ \mathbf{0} & L_b^{bq} \mathbf{1} + \mathbf{K}_b^{bq} & \mathbf{k}_b^{bp} & L_b^{bb} \mathbf{1} + \mathbf{K}_b^{bb} \end{bmatrix} \begin{bmatrix} T^{(q)} \\ \mathbf{v}^{(q)} \\ p^{(q)} \\ \lambda \end{bmatrix} \right) = \sum_{e=1}^{N_e} \left( \begin{bmatrix} f_T^{\text{ext.}(p)} \\ f_v^{\text{ext.}(p)} \\ 0 \\ f_b^{\text{ext}} \end{bmatrix} \right) \quad (47)
$$

where terms or submatrices $\square_2^{\gamma}$, $\square = \{M, L, \mathbf{K}, \mathbf{k}\}$, are defined relative to equations $z = \{T, v, p, b\}$, with shape functions associated to the couple of nodes $(x, y) = (p, q), (p, b), (b, q)$ or $(b, b)$ ($1 \le p, q \le 4$).

### 3.3.1. Linearization

We consider a fully implicit temporal discretization of the system (39). The residue at time $t + \Delta t$ is written as:

$$ \mathbf{r} = \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{conv}} - \mathbf{f}^{\text{int}} - \mathbf{M}\dot{\mathbf{q}} = \mathbf{0}. \quad (48) $$

The residue is linearized with a Newton-Raphson method at iteration $k$ of the computation leading to the following system:

$$ \mathbf{M}^{(k)} \dot{\mathbf{q}}^{(k+1)} + \mathbf{K}^{(k)} \delta \mathbf{q}^{(k)} = \mathbf{R}^{(k)} \quad (49) $$

where $\delta \mathbf{q}^{(k)} = \mathbf{q}^{(k+1)} - \mathbf{q}^{(k)}$ denotes the increment of the vector of degrees of freedom between two iterations and $\mathbf{R}^{(k)}$ a residue without the inertia terms $\mathbf{M}\dot{\mathbf{q}}$. The stiffness matrix $\mathbf{K}^{(k)}$ is defined as:

$$ \mathbf{K}^{(k)} = -\left. \frac{\partial \mathbf{R}}{\partial \mathbf{q}} \right|^{(k)} \quad (50) $$

and is assembled from contributions of submatrices pertaining to temperature, velocity and pressure degrees of freedom of vertex nodes and to degrees of freedom of the bubble node:

$$ \mathbf{K}_{xy} = -\frac{\partial \mathbf{R}_x}{\partial \mathbf{y}}, \quad (\mathbf{x}, \mathbf{y}) = (\mathbf{T}, \mathbf{v}, \mathbf{p}, \mathbf{\lambda}). \quad (51) $$

The system (49) is of the form:

$$
\begin{bmatrix} \mathbf{C} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{M}_{vv} & \mathbf{0} & \mathbf{M}_{vb} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{M}_{vb}^T & \mathbf{0} & \mathbf{M}_{bb} \end{bmatrix}^{(k)} \begin{Bmatrix} \dot{\mathbf{T}} \\ \dot{\mathbf{v}} \\ \dot{\mathbf{p}} \\ \dot{\mathbf{\lambda}} \end{Bmatrix}^{(k+1)} + \begin{bmatrix} \mathbf{K}_{TT} & \mathbf{K}_{Tv} & \mathbf{0} & \mathbf{K}_{Tb} \\ \mathbf{K}_{vT} & \mathbf{K}_{vv} & \mathbf{K}_{vp} & \mathbf{K}_{vb} \\ \mathbf{0} & \mathbf{K}_{pv} & \mathbf{0} & \mathbf{K}_{pb} \\ \mathbf{K}_{bT} & \mathbf{K}_{bv} & \mathbf{K}_{bp} & \mathbf{K}_{bb} \end{bmatrix}^{(k)} \begin{Bmatrix} \delta \mathbf{T} \\ \delta \mathbf{v} \\ \delta \mathbf{p} \\ \delta \mathbf{\lambda} \end{Bmatrix}^{(k)} = \begin{Bmatrix} \mathbf{R}_T \\ \mathbf{R}_v \\ \mathbf{R}_p \\ \mathbf{R}_b \end{Bmatrix}^{(k)}. \quad (52)
$$

The degrees of freedom of the bubble node $\mathbf{\lambda}$ are usually eliminated to save some computation cost by reducing the size of the system to be solved [9]. However, the presence of rates of bubble degrees of freedom $\mathbf{\lambda}$ makes this operation more complex. We therefore introduce below an approximation that allows to drop these rates.

### 3.3.2. Introduction of an extra approximation

Let us express the acceleration of the bubble node:

$$ \dot{\mathbf{v}}(\mathbf{x}^{(b)}, t) = \sum_{p=1}^{4} N^{(p)}(\mathbf{x}^{(b)}) \dot{\mathbf{v}}^{(p)}(t) + N^{(b)}(\mathbf{x}^{(b)}) \dot{\mathbf{\lambda}}(t). \quad (53) $$

The relation $\sum_{p=1}^4 N^{(p)} = 1$ and the fact that the bubble node is located at the centroid of the tetrahedron enable us to simplify expression (53). Indeed these elements imply that $N^{(p)}(\mathbf{x}^{(b)}) = \frac{1}{4}$. Moreover, the shape function associated to the bubble node is normalized so that $N^{(b)}(\mathbf{x}^{(b)}) = 1$. The acceleration at the bubble node can thus be rewritten as:

$$ \dot{\mathbf{v}}(\mathbf{x}^{(b)}, t) = \frac{1}{4} \sum_{p=1}^{4} \dot{\mathbf{v}}^{(p)}(t) + \dot{\mathbf{\lambda}}(t). \quad (54) $$

We shall now introduce an approximation. We assume that the bubble acceleration is roughly equal to the average of those of vertex nodes, which leads to:

$$ \dot{\mathbf{v}}(\mathbf{x}^{(b)}, t) \approx \frac{1}{4} \sum_{p=1}^{4} \dot{\mathbf{v}}^{(p)}(t). \quad (55) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Mesh of the viscometer.</figcaption></figure>

This assumption is valid when the acceleration field varies little on the element, and is compatible with a rigid body motion as detailed in the Appendix. Combination of expressions (54) and (55) leads to:

$$ \dot{\lambda} \approx 0. \qquad (56) $$

The elimination of the quantity $\lambda$ in Eq. (52) becomes then straightforward.

### 3.3.3. Solution

The elimination is performed at the element level accounting for (56), and time discretization is performed leading to the following linear system at each iteration:

$$ \mathbf{A}^{(k)} \delta \bar{\mathbf{q}}^{(k)} = \mathbf{B}^{(k)} \quad \text{with} \quad \left\{ \begin{array}{l} \mathbf{A}^{(k)} = \left( \frac{\bar{\mathbf{M}}^{(k)}}{\Delta t} + \bar{\mathbf{K}}^{(k)} \right) \\ \mathbf{B}^{(k)} = \bar{\mathbf{R}}^{(k)} - \bar{\mathbf{M}}^{(k)} \dot{\bar{\mathbf{q}}}^{(k)} \end{array} \right. \qquad (57) $$

where $\bar{\mathbf{q}}$ denotes the reduced vector of degrees of freedom of the system, defined as:

$$ \bar{\mathbf{q}}^T = \{\mathbf{T} \ \mathbf{v} \ \mathbf{p}\}. \qquad (58) $$

## 4. Comparison of analytical and numerical solutions

In this section, we present a comparison between analytical and numerical solutions for both cases of the viscometer problem considered in Section 2. The fluid domain is discretized with P1+/P1 finite elements, the mesh used is shown in Fig. 3. The inner and outer radii of the mesh are respectively set to $a = 0.1$ m and $b = 1$ m.

### 4.1. Unsteady mechanical state

In Eq. (13) for the unsteady mechanical state, the parameter $\lambda$ is arbitrary and arises from the change of variable (6) made in order to identify the differential equation (7) found to be a Bessel equation. A comparison between analytical and numerical solutions is made for two different values of this parameter to assess its influence.

The initial conditions are prescribed on each node, and velocities of the inner cylinder nodes are imposed according to a decreasing exponential. Thermal degrees of freedom have been fixed to zero at all nodes of the mesh. The numerical parameter values taken are the following: the rotation speed $\Omega_0$ (see Eq. (13)) of the inner radius is set to $\Omega_0 = 100 \text{ rad s}^{-1}$ and the kinematic viscosity is fixed to $\nu = \frac{\mu}{\rho} = 1 \text{ m}^2 \text{ s}^{-1}$. Actually, this high kinematic viscosity value does not correspond to any particular fluid; it is chosen in order to test more severely the dynamic behaviour of the formulation of the P1+/P1 finite element.

Fig. 4 presents the first plot obtained setting $\lambda = 1 \text{ m}^{-1}$, at various instants. The velocity field is plotted on a radial line. Points identified on the graph correspond to finite element solutions extracted from the nodes of the mesh on a radial line apparent in Fig. 3. We can observe a good correlation between analytical and numerical solutions.

It is worth noting that the case $\lambda = 1 \text{ m}^{-1}$ enables a physical interpretation. Indeed, this case is similar to a brake test of the inner cylinder, its rotation speed decreasing exponentially.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Transient velocity field for $\lambda = 1$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Transient velocity field for $\lambda = 5$.</figcaption></figure>

In Fig. 5, a plot obtained by setting $\lambda = 5 \text{ m}^{-1}$ is presented at various instants. Some comments are in order:

* As the arbitrary parameter $\lambda$ appears squared in the exponential of the solution (13), the transient velocity field decreases much more rapidly.

* The velocity field vanishes at some point. As $\lambda$ continues to grow, more and more oscillations of the velocity field will appear. This phenomenon arises from the fact that $\lambda$ appears in the argument of Bessel's functions in Eq. (13), which are oscillating functions.

* Again, we observe a good correlation between analytical and numerical solutions.

## 4.2. Unsteady thermal state-steady mechanical state

For this test case, the rotation speed of the inner cylinder is kept constant at $\Omega_0 = 100 \text{ rad s}^{-1}$. The thermal conductivity $k$, the density $\rho$ and the heat capacity $C$ are arbitrarily set so that the thermal diffusivity (Eq. (24)) be sufficiently high to test the thermal transient behaviour of the P1+/P1 finite element. Therefore, material parameters are fixed at $k = 200 \text{ W m}^{-1} \text{ K}^{-1}$, $\rho = 10 \text{ kg m}^{-3}$ and $C = 90 \text{ J kg}^{-1} \text{ K}^{-1}$.

Initial thermal and mechanical conditions extracted from the analytical solutions are prescribed on each node. The velocity field prescribed is given by Eq. (17), and the temperature field is given by Eq. (34) accounting for an initial condition such that $K$ be unity. The parameter $\lambda_i$ arising in (34) is a solution of Eq. (31).

Fig. 6 shows a comparison between analytical and numerical solutions for the unsteady thermal case at various instants, taking the first root $\lambda_1$. The temperature field is plotted on a radial line of the viscometer until its steady state. A good agreement between both solutions for the temperature field can be observed. The fully coupled thermofluid problem is nonlinear and may require a few iterations to converge (less than 10), thus the computation may last a little longer than those of the unsteady mechanical case, but never more than a few minutes.
<figure><img src="image_7.png" /><figcaption>Fig. 6. Transient temperature field with $λ_1 = 3.31394$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Transient temperature field with $λ_2 = 6.85758$.</figcaption></figure>

In Fig. 7, a comparison is presented taking the second root $λ_2$ of Eq. (31). This root arising in the argument of Bessel's functions in Eq. (34), the temperature field oscillates over the viscometer radius. In the same manner as for the unsteady mechanical case, in which the parameter $λ$ arises in (13), more and more oscillations of the temperature field will appear as a higher root $λ_i$ is used. A good agreement between analytical and numerical solutions is also observed.

## 5. Conclusion

In present Part I, extensions of the classical Couette viscometer solution in fluid mechanics have been developed to account for inertia effects and strong coupling between thermal and mechanical effects, the purpose of which is to assess numerical methods embedding such features. The benchmark tests developed are applied to the verification of a new formulation of the P1+/P1 finite element for transient laminar flow of incompressible fluids.

First, two analytical solutions of the Couette viscometer have been developed respectively for an unsteady mechanical state, and for a coupling (unsteady thermal state)–(steady mechanical state), considering a Newtonian fluid. Solutions have been 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; these solutions are expressed with Bessel's functions of the first and second kinds. More physically meaningful initial and boundary conditions could easily be considered, at the expense, however of greater complexity of the solution.

Subsequently, the mixed and unsteady temperature/velocity/pressure formulation of the P1+/P1 fluid finite element has been presented. This formulation is relevant for some welding applications in which a strong interaction between thermal and mechanical effects drives the formation of the weld, as in Friction Stir Welding-type processes. This formulation is an extension of Feulvarch et al. [11,12]'s previous works to the unsteady case including convection terms. The time discretization is based on the implicit (backward) Euler finite difference scheme leading to an iterative resolution. In order
to save computational cost, the bubble degrees of freedom are eliminated at the element level before the global procedure of the solution. An extra approximation made on the acceleration field enables to make this elimination, preserving the consistency of the semidiscrete equations with respect to rigid body motions.

Finally, comparisons between analytical and numerical solutions have been performed. These comparisons show a good agreement of both solutions developed, and attests to the good behaviour of the dynamic part and the thermomechanical coupling of the formulation of the finite element. However, it is worth noting that since the kinematics of the Couette viscometer associated to the incompressibility condition lead to a velocity and a temperature gradient orthogonal to each other, analytical solutions are developed with vanishing convection terms and therefore cannot claim to test the robustness of the implementation of these specific terms.

Part II will be devoted to the extension of the Couette viscometer problem to solid-type nonlinear behaviours, in the case of a purely quasi-static mechanical problem then in the case of a thermomechanical strongly-coupled problem, both for small and large geometry changes. These solutions will permit to assess a new formulation of a solid P1+/P1 finite element developed in the framework of large geometry changes with a temperature/velocity/pressure formulation.

## Appendix. Compatibility of the extra approximation with a rigid body motion

The approximation made on the bubble field is compatible with rigid body motions, in the sense that the weak form expressed with such motions reduces exactly to the Newton's second law. Indeed, the bubble node being located at the centroid of the tetrahedron, we have:

$$ \mathbf{x}^{(b)} = \frac{1}{4} \sum_{p=1}^{4} \mathbf{x}^{(p)}. \qquad (A.1) $$

For all rigid body motions, the bubble node remains located at the centroid of the element, therefore:

$$ \frac{d^2 \mathbf{x}^{(b)}}{dt^2} = \frac{1}{4} \sum_{p=1}^{4} \frac{d^2 \mathbf{x}^{(p)}}{dt^2} \Leftrightarrow \dot{\mathbf{v}}(\mathbf{x}^{(b)}) = \frac{1}{4} \sum_{p=1}^{4} \dot{\mathbf{v}}^{(p)}. \qquad (A.2) $$

This means that rigid body motions exactly satisfy approximations (55), which therefore does not introduce any error for such motions.

## References

[1] M. Couette, Études sur le frottement des liquides, Ann. Chim. Phys. Sér. 6 21 (1890) 433-510.

[2] G.G. Stokes, On the theories of the internal friction of fluids in motion, and of the equilibrium and motion of elastic solids, Trans. Cambridge Soc. 8 (1845) 287-305.

[3] D. Bernardin, Theoretical study of some transient Couette flows of viscoelastic fluids in inertial devices, J. Non-Newton. Fluid Mech. 88 (1999) 1-30.

[4] M. Kamran, M. Imran, M. Athar, M.A. Imran, On the unsteady rotational flow of fractional Oldroyd-B fluid in cylindrical domains, Meccanica 47 (2012) 573-584.

[5] D.N. Arnold, F. Brezzi, M. Fortin, A stable finite element for the Stokes equations, Calcolo 21 (1984) 337-344.

[6] P.M. Gresho, R.L. Sani, Incompressible Flow and the Finite Element Method: Advection-Diffusion and Isothermal Laminar Flow, Wiley, 1998.

[7] J.L. Chenot, F. Bay, An overview of numerical modelling techniques, J. Mater. Process. Technol. 80-81 (1998) 8-15.

[8] M. Bellet, V.D. Fachinotti, ALE method for solidification modelling, Comput. Methods Appl. Mech. Engrg. 193 (2004) 4355-4381.

[9] M. Bellet, O. Jaouen, I. Poitrault, An ALE-FEM approach to the thermomechanics of solidification processes with application to the prediction of pipe shrinkage, Intl. J. Numer. Methods Heat Fluid Flow 15 (2005) 120-142.

[10] O. Basset, Simulation numérique d'écoulements multi-fluides sur grille de calcul. Ph.D. thesis, École des Mines de Paris, 2006 (in French).

[11] É. Feulvarch, Y. Gooroochurn, F. Boitout, J.M. Bergheau, 3D modelling of thermofluid flow in Friction Stir Welding, in: S.A. David, T. DebRoy, J.C. Lippold, H.B. Smartt, J.M. Vitek (Eds.), Proc. 7th Int. Conf. on Trends in Welding Research, Pine Mountain, Georgia, USA, 2005, pp. 261-266.

[12] É. Feulvarch, F. Boitout, J.M. Bergheau, Friction Stir Welding: modélisation de l'écoulement de la matière pendant la phase de soudage, Eur. J. Comp. Mech. 16 (2007) 865-887 (in French).

[13] É. Feulvarch, J.C. Roux, J.M. Bergheau, A simple and robust moving mesh technique for the finite element simulation of Friction Stir Welding, J. Comput. Appl. Math. 246 (2013) 269-277.

[14] SYSWELD® (2007). ESI Group. User's manual.

[15] T. Heuzé, J.B. Leblond, J.M. Bergheau, É. Feulvarch, A finite element for laminar flow of incompressible fluids with inertia effects and thermomechanical coupling, Eur. J. Comp. Mech. 19 (2010) 293-304.

[16] T. Heuzé, Modélisation des couplages fluide/solide dans les procédés d'assemblage à haute température, Ph.D. thesis, Université Pierre et Marie Curie, Paris, 2011 (in French).

[17] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Second ed., Clarendon press -cop., Oxford, 1953.

[18] A.N. Brooks, T.J.R. Hughes, Streamline UpWind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199-259.

[19] T. Belytschko, W.K. Liu, B. Moran, Nonlinear Finite Elements for Continua and Structures, Wiley, 2000.

[20] J.M. Bergheau, R. Fortunier, Finite Element Simulation of Heat Transfert, Wiley, 2008.