# Benchmark tests based on the Couette viscometer

Thomas Heuzé, Hussein Amin El Sayed, Jean-baptiste Leblond, Jean-michel Bergheau

► To cite this version:

Thomas Heuzé, Hussein Amin El Sayed, Jean-baptiste Leblond, Jean-michel Bergheau. Benchmark tests based on the Couette viscometer: II: Thermo-elasto-plastic solid behaviour in small and large strains. Computers & Mathematics with Applications, 2014, 67 (8), pp.1482-1496. 10.1016/j.camwa.2014.02.010. hal-01667115

HAL Id: hal-01667115

https://hal.science/hal-01667115

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—II: Thermo-elasto-plastic solid behaviour in small and large strains

Thomas Heuzé<sup>a,*</sup>, Hussein Amin-El-Sayed<sup>b</sup>, Jean-Baptiste Leblond<sup>c</sup>, Jean-Michel Bergheau<sup>b</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 and to use the extended solutions to assess new finite elements. Part I was devoted to the case of laminar flow of incompressible fluids with inertia effects and thermomechanical coupling. The present Part II focuses on solid-type nonlinear behaviours; we address the cases of elastic-plastic and thermo-elastic-plastic von Mises materials, both in small and large strains. The extended solutions permit to assess a new formulation of a mixed P1 + /P1 finite element in solid mechanics, in a temperature/velocity/pressure formulation coupled with an implicit (backward) Euler algorithm in time. The verification evidences a good behaviour of the solid finite element.

## 1. Introduction

Part I was devoted to the viscometer problem in the case of fluid-type behaviours; extended solutions were defined accounting for inertia effects and strong coupling between thermal and mechanical effects. Though the viscometer is generally considered as a typical fluid mechanics apparatus or problem, it seems natural to keep the same problem to develop benchmark tests for solid-type nonlinear behaviours. Extensions of the problem to solid materials are therefore the topic of the present Part II.

In Section 2, we consider the case of an elastic-plastic material obeying the von Mises criterion with isotropic hardening. The elastic-plastic problem is first studied in Section 2.1 in the small strain framework. It is shown that the displacement field is found by integrating a function defined by the elastic-plastic constitutive law. The load curve of the structure and the unloading stage are then studied in Sections 2.3 and 2.4. It is shown that for the Couette viscometer the unloading leads, somewhat surprisingly, to a vanishing residual stress field. An extension to the large strain framework is then presented in Section 2.5.

<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> Université de Lyon, ENISE, LTDS, UMR 5513 CNRS, 58 rue Jean Parot, 42023 Saint-Etienne, Cedex 2, France</footnote>

<footnote><sup>c</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>* Corresponding author.
E-mail address: thomas.heuze@ec-nantes.fr (T. Heuzé).</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Parameterization of the Couette viscometer problem.</figcaption></figure>

In Section 3, the Couette viscometer problem is addressed in thermo-elasto-plasticity. The solution permits to account for the generation of heat through plastic dissipation within the structure. The solution developed is valid only during the phase when the structure is entirely elastic-plastic.

The purpose of the extension of the Couette viscometer problem to solid behaviours is to assess a formulation in solid mechanics of the P1 + /P1 finite element initially introduced by Arnold et al. [1]. This element is developed in Section 4 in a new, fully coupled temperature/velocity/pressure formulation in the large strain framework, combined with an implicit (backward) Euler algorithm in time. The bubble node embedded in this element enables to ensure, with the mixed formulation, plastic incompressibility. The resolution is performed on the unknown configuration at time $t + \Delta t$. This element is one feature of the modelling of the Friction Stir Spot Welding process [2] recently implemented in the finite element code SYSWELD® [3].

The verification of the solid P1 + /P1 finite element is subsequently performed using the reference solutions developed, and presented in Section 5. First, the mechanical behaviour of the finite element is tested with respect to the solution developed in Section 2.5 in the large strain framework. The comparison shows a good agreement between numerical and analytical solutions. Second, the thermomechanical behaviour of the finite element is tested with respect to the solution developed in Section 3. Comparisons on thermal and mechanical fields also show a good agreement between numerical and analytical solutions.

## 2. The Couette viscometer problem for an elastic-plastic von Mises material with strain hardening

### 2.1. Small strain framework

We consider the Couette viscometer problem schematized in Fig. 1. The inner and outer cylinder radii are denoted $a$ and $b$ respectively. We shall consider in the sequel the inner cylinder as fixed and the outer cylinder as driven:

$$ \begin{aligned} u_{\theta}(r=a) &= 0 \\ u_{\theta}(r=b) &= u_{\theta}(b) \end{aligned} \qquad (1) $$

where $u_{\theta}$ denotes the orthoradial displacement. The aim of this choice, which differs from that made in Part I for fluid-type behaviours, is to warrant a positive shear stress $\sigma_{r\theta} \equiv \tau$.

We consider in this section an elastic-plastic material obeying a von Mises criterion with isotropic hardening. The yield stress is supposed not to saturate, that is to increase indefinitely with the plastic strain. The study is performed in the small strain framework with an increasing monotonic loading, starting from a natural initial state. Since an elastic-plastic constitutive law is used, three regimes occur during the loading:

* the fully elastic regime: the structure deforms elastically.

* the mixed elastic-plastic/elastic regime: an elastic-plastic crown appears adjacent to the inner radius, surrounded by an elastic crown.

* the fully elastic-plastic regime: the whole structure is elastic-plastic.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Elastic-plastic crown radius c in the mixed elastic-plastic/elastic regime.</figcaption></figure>

The stress field within the structure is of the form:

$$ \sigma = \tau(\mathbf{e}_r \otimes \mathbf{e}_\theta + \mathbf{e}_\theta \otimes \mathbf{e}_r) \qquad (2) $$

where $\sigma$ denotes the Cauchy stress tensor and $\tau$ a shear stress component. The tangential equilibrium equation:

$$ \frac{d\tau}{dr} + \frac{2\tau}{r} = 0 \qquad (3) $$

implies that this shear stress must be of the form:

$$ \tau = \frac{A}{r^2} \qquad (4) $$

where A is a constant which can be determined from appropriate boundary conditions. This constant may be rewritten as:

$$ A = kc^2 \qquad (5) $$

leading to the following expression of the shear stress:

$$ \tau = k \frac{c^2}{r^2} \qquad (6) $$

where k is the initial shear yield stress and c a parameter homogeneous to a length, which can serve as a loading parameter. Plasticity occurs at a given point r if $\tau \ge k$ or equivalently $c \ge r$. Thus, it follows that if:

* $c < a$, the structure is completely elastic,

* $a \le c < b$, the structure is in a mixed elastic-plastic/elastic regime,

* $c \ge b$, the structure is completely elastic-plastic.

In the second case, the constant c denotes the radius of the inner elastic-plastic crown, as shown in Fig. 2.

It remains to calculate the displacement field to complete the solution. This will be done by using the following equations:

$$ \tau = k \frac{c^2}{r^2} \quad (\text{equilibrium}) \\ \tau = f(\gamma) \quad (\text{constitutive law}) \qquad (7) $$

$$ \gamma = 2\epsilon_{r\theta} = \frac{du_{\theta}}{dr} - \frac{u_{\theta}}{r} \quad (\text{definition of strain}). $$

The function $f(\gamma)$ appearing in Eq. (7) is supposed to be strictly increasing.

The constitutive law (7) combined with the expression (7) of the stress permits to express the strain as a function of the load parameter c:

$$ \gamma = f^{-1}\left(k \frac{c^2}{r^2}\right) \equiv g\left(\frac{r}{c}\right) \qquad (8) $$
<figure><img src="image_4.png" /><figcaption>Fig. 3. Strain hardening law.</figcaption></figure>

where g is a function defined from the inverse of f. The definition (7) of the strain implies that:

$$ \gamma = \frac{du_{\theta}}{dr} - \frac{u_{\theta}}{r} = r \frac{d}{dr} \left( \frac{u_{\theta}}{r} \right) = g \left( \frac{r}{c} \right). \qquad (9) $$

Thus, the displacement field is obtained through integration:

$$ u_{\theta}(r) = r \int_{a}^{r} \frac{g(r'/c)}{r'} dr'. \qquad (10) $$

The displacement of the outer cylinder $u_{\theta}(b)$ is therefore given by:

$$ u_{\theta}(b) = b \int_{a}^{b} \frac{g(r'/c)}{r'} dr'. \qquad (11) $$

It is worth noting that $u_{\theta}(b)$ can be explicitly expressed as a function of c, whereas the converse is not true. This justifies the choice of c rather than $u_{\theta}(b)$ as a load parameter.

## 2.2. Specialization of the solution

The general solution for the displacement field is given by Eq. (10), provided that the function g is defined from an increasing function f. It is interesting to specialize this solution to a given constitutive law in order to derive an explicit expression of the displacement field, allowing for a comparison between numerical and analytical solutions. We shall consider a linear hardening law in the sequel, see Fig. 3.

Explicit expressions of the displacement field will be derived for each regime occurring during the monotonic loading of the structure.

### 2.2.1. Fully elastic regime ($c < a$)

In the linear regime, the function g is noted $g^e$ and is defined by:

$$ g^e \left( \frac{r}{c} \right) = \gamma = \frac{\tau(r)}{G} = \frac{k c^2}{G r^2} \qquad (12) $$

where G denotes the elastic shear modulus. This leads upon calculation of the integral of Eq. (10) to:

$$ u_{\theta}(r) = \frac{kr}{2G} \frac{c^2}{a^2} \left( 1 - \frac{a^2}{r^2} \right). \qquad (13) $$

In particular, the parameter $u_{\theta}(b)$ is related to c through the relation:

$$ u_{\theta}(b) = \frac{kb}{2G} \frac{c^2}{a^2} \left( 1 - \frac{a^2}{b^2} \right). \qquad (14) $$

### 2.2.2. Mixed elastic-plastic/elastic regime ($a \le c < b$)

In the mixed elastic-plastic/elastic regime, an elastic-plastic crown is surrounded by an elastic crown, thus the calculation of the integral (10) is performed separately in these two zones.
### Elastic-plastic zone

Since the initial state is supposed to be stress-free and the loading monotonically increasing, the elasticity law, the von Mises criterion and the expression of the cumulated plastic strain $\varepsilon_{eq}$ read:

$$ 
\begin{cases}
\tau = G(\gamma - \gamma^p) \\
\tau = k + h \frac{\varepsilon_{eq}}{\sqrt{3}} \\
\varepsilon_{eq} = \frac{\gamma^p}{\sqrt{3}}
\end{cases}
\quad (15) $$

where h denotes the strain hardening parameter in simple tension. Elimination of $\gamma^p$ and $\varepsilon_{eq}$ in Eqs. (15) gives:

$$ \gamma = \frac{h + 3G}{hG} \tau - \frac{3k}{h}. \quad (16) $$

Introducing Eq. (16) into (10), one gets the displacement field within the elastic-plastic zone:

$$ u_{\theta}(r) = kr \left[ \frac{h + 3G}{hG} \frac{c^2(r^2 - a^2)}{2a^2r^2} - \frac{3}{h} \ln\left(\frac{r}{a}\right) \right]. \quad (17) $$

The cumulated plastic strain is in turn obtained by combining Eqs. (15) 2 and (6):

$$ \varepsilon_{eq}(r) = \frac{\sqrt{3}k}{h} \left( \frac{c^2}{r^2} - 1 \right). \quad (18) $$

### Elastic zone

The elastic crown surrounds the elastic-plastic one, therefore the integral in Eq. (10) must be split in two to account for both the elastic-plastic and elastic zones:

$$ u_{\theta}(r) = r \left[ \int_{a}^{c} \frac{g^p(r'/c)}{r'} dr' + \int_{c}^{r} \frac{g^e(r'/c)}{r'} dr' \right] \quad (19) $$

where g takes the following expressions, noted $g^e$ and $g^p$ in the elastic and elastic-plastic zones respectively:

$$ * g^e(r'/c) = \frac{k c^2}{G r'^2} \quad \text{for } c \le r' \le r, \quad (20) $$

$$ * g^p(r'/c) = \frac{h + 3G}{hG} \tau(r') - \frac{3k}{h} \quad \text{for } a \le r' \le c. \quad (21) $$

The calculation of the integrals provides the following solution:

$$ u_{\theta}(r) = kr \left[ \frac{h + 3G}{hG} \frac{c^2 - a^2}{2a^2} - \frac{3}{h} \ln\left(\frac{c}{a}\right) + \frac{1}{2G} \left(1 - \frac{c^2}{r^2}\right) \right]. \quad (22) $$

In particular, the parameter $u_{\theta}(b)$ can be expressed as a function of the plastic crown radius c:

$$ u_{\theta}(b) = kb \left[ \frac{h + 3G}{hG} \frac{c^2 - a^2}{2a^2} - \frac{3}{h} \ln\left(\frac{c}{a}\right) + \frac{1}{2G} \left(1 - \frac{c^2}{b^2}\right) \right]. \quad (23) $$

### 2.2.3. Fully elastic-plastic regime ($c \ge b$)

When the whole structure is elastic-plastic, the sole function $g^p$ appears in the integral (10), thus the expression of the displacement field is given by (17). In particular, the parameter $u_{\theta}(b)$ is connected to the radius c by the relation:

$$ u_{\theta}(b) = kb \left[ \frac{h + 3G}{hG} \frac{c^2}{2a^2} \left( 1 - \frac{a^2}{b^2} \right) - \frac{3}{h} \ln \left( \frac{b}{a} \right) \right]. \quad (24) $$

The cumulated plastic strain is still given by Eq. (18).
<figure><img src="image_5.png" /><figcaption>Fig. 4. Load-displacement curve of the structure.</figcaption></figure>

## 2.3. Load-displacement curve of the structure

Dual force and displacement parameters $Q, q$ can be defined in such a way that the power of external forces read:

$$ \mathcal{P}_{\text{ext}} = \mathbf{Q} \cdot \dot{\mathbf{q}}. \qquad (25) $$

In the case of the Couette viscometer with the loading considered, these parameters may be identified as:

$$ \begin{aligned} Q &= M_z(b) \\ q &= \frac{u_{\theta}(b)}{b} \end{aligned} \qquad (26) $$

where the torque per unit thickness $M_z(b)$ is defined as:

$$ M_z(b) = \mathbf{M}(O) \cdot \mathbf{z} = \left[ \int b \mathbf{e}_r \times (\boldsymbol{\sigma} \cdot \mathbf{e}_r) b d\theta \right] \cdot \mathbf{z} = 2\pi b^2 \tau(b). \qquad (27) $$

This torque can be calculated for each regime occurring during the loading:

* **fully elastic regime**
Combining formulae (12)–(14) and (27), one gets:

$$ M_z(b) = \frac{4\pi G a^2 b^2 u_{\theta}(b)}{b^2 - a^2} \qquad (28) $$

* **mixed elastic-plastic/elastic regime**
One gets similarly:

$$ M_z(b) = 2\pi k c^2 \qquad (29) $$

where c is a function of $u_{\theta}(b)$ defined implicitly by relation (23).

* **fully elastic-plastic regime**
The plastic crown radius c is given here by (24), thus the expression of the torque reads:

$$ M_z(b) = \frac{4\pi k a^2 b^2}{b^2 - a^2} \frac{hG}{h + 3G} \left[ \frac{u_{\theta}(b)}{kb} + \frac{3}{h} \ln\left(\frac{b}{a}\right) \right]. \qquad (30) $$

The load-displacement curve of the structure is shown in Fig. 4. The structure remains completely elastic up to the load of incipient plasticity $M_z^{(1)}(b)$. Subsequently, a plastic crown gradually develops from the inner cylinder and extends radially until the whole structure becomes elastic-plastic, for the second torque $M_z^{(2)}(b)$. Afterwards, the torque increases linearly with the kinematic parameter $u_{\theta}(b)$, as a result of strain hardening.
## 2.4. Unloading

It is interesting to study an unloading of the structure from some completely elastic-plastic state and calculate residual stresses. The unloading leads to a zero force parameter $Q$. In the case of the Couette viscometer, the torque on the outer cylinder thus goes to zero, and so does also the shear stress at $r = b$. But it has been shown using the equilibrium equation (7)$_1$ that the shear stress varies as $1/r^2$. Therefore, the residual shear stress vanishes over the whole structure when it is unloaded.

This property is of course tied to the special nature of the problem considered and by no means general. More precisely, the plastic strain field is always of the form $\epsilon_{r\theta}^p = \gamma^p(r)/2$, other $\epsilon_{ij}^p = 0$. Such a strain field is always compatible, because the differential equation $\frac{du_{\theta}^p}{dr} - \frac{u_{\theta}^p}{r} = \gamma^p(r)$ (where $u_{\theta}^p$ is some unknown “plastic displacement”) always admits a solution; it is a well-known property of plasticity theory that the residual stresses after unloading are then necessarily zero [4,5].

## 2.5. Extension to the large strain framework

The Couette viscometer problem for an elastic-plastic material has been considered in the small strain assumption up to now. Its extension to the large strain framework is important for the assessment of a new formulation of the P1 + /P1 finite element presented in Section 4.

The plasticity equations in large strain are written here in the Eulerian form; the linearized strain rate tensor $\dot{\epsilon}$ is replaced by the Eulerian strain rate tensor $\mathbf{D}$. The latter tensor is decomposed additively into elastic and plastic parts $\mathbf{D}^e, \mathbf{D}^p$. The elastic strain rate tensor $\mathbf{D}^e$ is linked to an objective stress rate through some hypoelastic law, while the plastic strain rate tensor $\mathbf{D}^p$ is given by the flow rule associated to von Mises's criterion.

The presence of some objective derivative of the stress tensor in the hypoelastic law makes the solution more complex. However, if the additional terms associated to this objective derivative are neglected, the solution developed for small strains remains unchanged provided the tangential displacement $u_{\theta}$ is replaced by the curved arc length swept $l$. This point is detailed in the Appendix.

# 3. The Couette viscometer problem in thermo-elasto-plasticity

We develop in this section a solution of the Couette viscometer problem in thermo-elasto-plasticity, accounting for the generation of heat through plastic dissipation within the structure. This solution allows for comparisons with the P1 + /P1 solid finite element developed with a strong thermomechanical coupling presented in Section 4. The mechanical part of the solution corresponds to an elastic-plastic constitutive law with the von Mises criterion and a linear isotropic strain hardening (Fig. 3), already considered in both small and large strain in Section 2. The thermal part of the solution is governed by the equation of energy conservation:

$$ \rho C \frac{\partial T}{\partial t} = \nabla \cdot (\lambda \nabla T) + S \quad \forall \mathbf{x} \in \Omega \quad (31) $$

supplemented with appropriate initial and boundary conditions. The parameters $\rho$, $C$ and $\lambda$ here denote the density, the heat capacity and the thermal conductivity, respectively, and $S$ is a source term. The kinematics of the Couette viscometer permit to simplify Eq. (31) into:

$$ \rho C \frac{\partial T}{\partial t} = \frac{\lambda}{r} \frac{\partial}{\partial r} \left( r \frac{\partial T}{\partial r} \right) + S. \quad (32) $$

The strong thermomechanical coupling consists on the one hand of the effects of the mechanical dissipation appearing in Eqs. (31) and (32) in the source term $S$, and on the other hand of thermal dilation and temperature dependent mechanical parameters. To facilitate the solution, the influence of the thermal part of the solution upon its mechanical part is not accounted for here, that is dilation effects are neglected and mechanical material parameters are fixed, independently of temperature.

With these assumptions, the mechanical dissipation is given by:

$$ S = \beta \tau \dot{\gamma}^p \quad (33) $$

where $\beta$ denotes the Quinney-Taylor coefficient accounting for the fact that only a portion of the plastic power is dissipated into heat; this coefficient is usually set to 0.9. The mechanical dissipation $S$ is here time-dependent, thus the thermal part of the solution depends on its mechanical part determined a priori. Assuming the structure to be entirely elastic-plastic, the shear stress $\tau$ is given by Eq. (6) and the plastic rate of slip $\dot{\gamma}^p$ by the expression $3\dot{\tau}/h$ resulting from Eqs. (15)$_2$ and (15)$_3$; the mechanical dissipation can then be expressed as:

$$ S = \frac{6\beta k^2 c^3 \dot{c}}{h r^4}. \quad (34) $$
To solve the partial differential equation (32), it is convenient to assume that both the plastic crown radius c and the temperature vary exponentially in time, in order to allow for a separation of space and time coordinates. Let the evolution of the plastic crown radius be thus given by:

$$c(t) = c_0 \exp(\alpha t) \tag{35}$$

where $\alpha$ is a constant homogeneous to the inverse of a time. The mechanical dissipation is then expressed as:

$$S = \frac{6\beta k^2 c_0^4 \alpha}{h} \frac{\exp(4\alpha t)}{r^4} \tag{36}$$

The temperature is now assumed to be of the form:

$$T(r, t) = f(r) \exp(4\alpha t) \tag{37}$$

The partial differential equation (32) reduces to the following inhomogeneous ordinary differential equation on the unknown function $f(r)$:

$$f''(r) + \frac{f'(r)}{r} - \mu^2 f(r) = - \frac{6\beta k^2 c_0^4 \alpha}{h\lambda} \frac{1}{r^4} \quad \text{where } \mu^2 = \frac{4\rho C\alpha}{\lambda} \tag{38}$$

To solve Eq. (38), we first perform a change of variable and function:

$$x = \mu r, \quad F(x) = f(r) \tag{39}$$

where since $x$ is dimensionless, $\mu$ is homogeneous to the inverse of a length. Eq. (38) can be rewritten, accounting for (39), as:

$$F''(x) + \frac{F'(x)}{x} - F(x) = -\frac{B}{x^4} \tag{40}$$

where $B$ is a constant defined by:

$$B = \frac{6\beta k^2 c_0^4 \alpha \mu^2}{h\lambda} \tag{41}$$

Eq. (40) can be identified to a modified inhomogeneous Bessel's equation [6] at order zero.

We look for a solution of Eq. (40) in the form:

$$F(x) = I_0(x)G(x) \tag{42}$$

where $I_0(x)$ is the modified Bessel function of the first kind at order zero, and $G(x)$ an unknown function. Introducing (42) into (40), and noting that $I_0(x)$ satisfies the homogeneous equation associated to (40), one gets:

$$G''(x) + \left( \frac{1}{x} + \frac{2I_0'(x)}{I_0(x)} \right) G'(x) = - \frac{B}{I_0(x)x^4} \tag{43}$$

By setting $H(x) = G'(x)$, Eq. (40) is reduced to the first order differential equation

$$H'(x) + \left( \frac{1}{x} + \frac{2I_0'(x)}{I_0(x)} \right) H(x) = - \frac{B}{I_0(x)x^4} \tag{44}$$

The solution of the homogeneous equation associated to (44) is obtained through separation of variables and is given by:

$$H(x) = \frac{K}{xI_0^2(x)} \tag{45}$$

where $K$ is a constant. The method of variation of the constant then permits to determine the solutions of Eq. (44); one thus finds that:

$$H(x) = \frac{B\phi(x) - C_1}{xI_0^2(x)} \tag{46}$$

where $C_1$ is a constant and the function $\phi(x)$ is defined by:

$$\phi(x) = \int_x^s \frac{I_0(u)}{u^3} du \quad \forall x > 0. \tag{47}$$
The upper bound s of the integral defining $\phi(x)$ may be set arbitrarily (changing it just changes the value of the constant $C_1$). However, since the integrand of $\phi(x)$ goes to infinity for $u \to +\infty$, s should be ascribed some finite value.

To calculate the function $G(x)$, the expression (46) of the function $H(x)$ needs to be integrated. Use will be made here of the following identity which can be established by combining formulas 9.6.15 and 9.6.27 of Abramowitz and Stegun [7]:

$$ \left(\frac{K_0}{I_0}\right)'(x) = -\frac{1}{xI_0^2(x)} \qquad (48) $$

where $K_0(x)$ is the modified Bessel function of the second kind at order zero. According to Eq. (48), the expression of $H(x)$ may be written in the form:

$$ H(x) = (-B\phi(x) + C_1) \left(\frac{K_0}{I_0}\right)'(x). \qquad (49) $$

The function $G(x)$ is then obtained through integration. The second term of Eq. (49) is directly integrated whereas the first one is integrated by parts. The solution of Eq. (40) finally reads:

$$ F(x) = (C_2 + B\psi(x))I_0(x) + (C_1 - B\phi(x))K_0(x) \qquad (50) $$

where $C_2$ is an additional constant and the function $\psi(x)$ is defined by:

$$ \psi(x) = \int_{x}^{+\infty} \frac{K_0(u)}{u^3} du \quad \forall x > 0. \qquad (51) $$

Taking an infinite upper bound here does not raise any problem since $K_0(u)$ goes to zero when $u$ goes to infinity.

Boundary conditions are defined in such a way as to lead to the simplest possible expression of the solution field. We thus prescribe homogeneous Dirichlet conditions on both cylinders:

$$ \begin{cases} T(r=a, t)=0 \Rightarrow f(a)=0 \Rightarrow F(\mu a)=0 \\ T(r=b, t)=0 \Rightarrow f(b)=0 \Rightarrow F(\mu b)=0. \end{cases} \qquad (52) $$

Note that other couples of boundary conditions, such as zero fluxes imposed on the cylinders, could also be prescribed. Solving the system of Eqs. (52) leads to the following values of the constants $C_1$ and $C_2$:

$$ C_1 = \frac{B[-I_0(\mu a)K_0(\mu b)\phi(\mu b) + I_0(\mu b)(K_0(\mu a)\phi(\mu a) - I_0(\mu a)(\psi(\mu a) - \psi(\mu b)))]}{I_0(\mu b)K_0(\mu a) - I_0(\mu a)K_0(\mu b)} \qquad (53) $$

$$ C_2 = \frac{B[-I_0(\mu a)K_0(\mu b)\psi(\mu a) + K_0(\mu a)(K_0(\mu b)(\phi(\mu a) - \phi(\mu b)) + I_0(\mu b)\psi(\mu b))]}{I_0(\mu a)K_0(\mu b) - I_0(\mu b)K_0(\mu a)}. $$

Therefore the transient temperature field in the coupled thermo-elasto-plastic viscometer problem is given by:

$$ \begin{align} T(r,t) &= [(C_2 + B\psi(\mu r))I_0(\mu r) + (C_1 - B\phi(\mu r))K_0(\mu r)] \exp(4\alpha t) \nonumber \\ &\quad \text{with } \alpha = \frac{\lambda\mu^2}{4\rho C} \quad (\mu \text{ arbitrary parameter}) \nonumber \\ B &= \frac{6\beta k^2 c_0^4 \alpha \mu^2}{h\lambda}, \quad C_1, C_2 \text{ given by (53), } \phi(x) \text{ by (47), } \psi(x) \text{ by (51)} \tag{54} \\ T(r=a, t) &= 0 \nonumber \\ T(r=b, t) &= 0 \nonumber \\ T(r,t=0) &= T_0(r). \nonumber \end{align} $$

# 4. Mixed temperature/velocity/pressure formulation for the P1 + /P1 solid finite element

## 4.1. Weak formulation

In welding applications, geometrical nonlinearities such as large rotations often occur, and govern the deformed shape of the structure [8–16]. Moreover, certain welding processes such as Friction Stir Welding also involve large strains [17,18]. In addition, Friction Stir Welding involves a strong coupling between thermal and mechanical effects. Indeed, the heat generated by friction and the motion of the material generates a mix which joins the parts to be welded after cooling.

We propose a special formulation in solid mechanics within the large displacement/strain framework (Eulerian setting) which includes a strong thermomechanical coupling. 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 (55) $$

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

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

In a mixed formulation, the weak form of the problem is as follows: given the body forces $\mathbf{f}$, the heat source $r$, the prescribed tractions $\mathbf{F}_d$ on $\partial \Omega_F$, and the prescribed 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 } (T, \mathbf{v}, p) \in (\mathcal{T}_{\text{ad}} \times \mathcal{V}_{\text{ad}} \times \mathcal{P}_{\text{ad}}), \quad \forall t \in [0, T], \text{ such that } \forall (T^*, \mathbf{v}^*, p^*) \in (\mathcal{T}_{\text{ad}}^0 \times \mathcal{V}_{\text{ad}}^0 \times \mathcal{P}_{\text{ad}}), \\ \quad \begin{aligned} & -\int_{\Omega} k \nabla T \cdot \nabla T^* d\Omega + \int_{\partial \Omega_q} \phi^d T^* dS + \int_{\Omega} \beta \boldsymbol{\sigma} : \mathbf{D}^p T^* d\Omega + \int_{\Omega} r T^* d\Omega = \int_{\Omega} \rho C \frac{\partial T}{\partial t} T^* d\Omega \\ & -\int_{\Omega} \boldsymbol{\sigma} : \mathbf{D}^* d\Omega + \int_{\partial \Omega_F} \mathbf{F}_d \cdot \mathbf{v}^* dS + \int_{\Omega} \mathbf{f} \cdot \mathbf{v}^* d\Omega = 0 \end{aligned} \\ \quad \int_{\Omega} p^* \text{tr} \mathbf{D}^p d\Omega = 0 \\ T(\mathbf{x}, t=0) = T_0(\mathbf{x}) \\ \mathbf{v}(\mathbf{x}, t=0) = \mathbf{v}_0(\mathbf{x}) \end{array} \right. \qquad (57) $$

where $\beta$ denotes the Quinney-Taylor coefficient introduced in Section 3, and $\mathbf{D}^p$ the plastic part of the strain rate. The functional spaces involved here are defined as

$$ \begin{align} \mathcal{V}_{\text{ad}} &= \{\mathbf{v}(\mathbf{x}, t) \in [H^1(\Omega)]^3 | \mathbf{v}(\mathbf{x}, t) = \mathbf{v}^d(\mathbf{x}, t) \text{ on } \partial \Omega_v\} \notag \\ \mathcal{V}_{\text{ad}}^0 &= \{\mathbf{v}^* \in [H^1(\Omega)]^3 | \mathbf{v}^* = \mathbf{0} \text{ on } \partial \Omega_v\} \notag \\ \mathcal{T}_{\text{ad}} &= \{T(\mathbf{x}, t) \in H^1(\Omega) | T(\mathbf{x}, t) = T^d(\mathbf{x}, t) \text{ on } \partial \Omega_{\theta}\} \tag{58} \\ \mathcal{T}_{\text{ad}}^0 &= \{T^* \in H^1(\Omega) | T^* = 0 \text{ on } \partial \Omega_{\theta}\} \notag \\ \mathcal{P}_{\text{ad}} &= \{p^* \in L^2(\Omega)\} \notag \end{align} $$

where $\mathbf{v}^d(\mathbf{x}, t)$ and $T^d(\mathbf{x}, t)$ are the velocity and temperature prescribed on $\partial \Omega_v$ and $\partial \Omega_{\theta}$. The mechanical inertia terms have been neglected in Eqs. (57) since the deformation process may safely be considered as quasi-static.

The Cauchy stress tensor is decomposed into its spherical and traceless parts:

$$ \boldsymbol{\sigma} = \mathbf{p}\mathbf{1} + \mathbf{s} \qquad (59) $$

where $\mathbf{p}$ is the opposite of the hydrostatic pressure, **1** the second-rank identity tensor, and **s** the stress deviator. For solid-type behaviours, the hydrostatic pressure is not determined by solving the initial boundary value problem like in fluid mechanics, but is related to the trace of the elastic part of the strain rate $\mathbf{D}^e$ through the elastic law:

$$ \mathrm{tr} \mathbf{D}^e = \frac{3(1 - 2\nu)}{E} \dot{\mathbf{p}} = \frac{\dot{\mathbf{p}}}{\kappa} \qquad (60) $$

where $\kappa$ denotes the elastic bulk modulus. Using the decomposition of the strain rate and relation (60), the expression of the plastic incompressibility condition, written in weak form in (57), becomes:

$$ \int_{\Omega} p^* \left( \nabla \cdot \mathbf{v} - \frac{\dot{\mathbf{p}}}{\kappa} - \text{tr} \dot{\boldsymbol{\varepsilon}}^{\text{th}} \right) d\Omega = 0, \quad \forall p^* \in \mathcal{P}_{\text{ad}} \qquad (61) $$

where $\boldsymbol{\varepsilon}^{\text{th}}$ denotes the thermal strain tensor.

## 4.2. Finite element formulation

### 4.2.1. Semidiscrete equations

The weak form (57) is discretized with the P1 + /P1 finite element. Basic features of this finite element have already been presented in Part I. The finite element discretization leads to the following system of semidiscrete equations:

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

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

$$ \mathbf{q}^T = \{\mathbf{T} \mathbf{v} \mathbf{p} \boldsymbol{\lambda}\}. \qquad (63) $$
The mass matrix **M** consists only of the heat capacity matrix, the mechanical inertia terms being neglected. The internal and external forces **f**<sup>int</sup>, **f**<sup>ext</sup> are defined by:

$$
\mathbf{f}^{\text{int}} = \sum_{e=1}^{N_e} \left( \begin{array}{l} \left\{ \begin{aligned} f_T^{\text{int},(p)} &= -\int_{\Omega_e} \nabla N^{(p)} \cdot \mathbf{hd}\Omega \\ \mathbf{f}_v^{\text{int},(p)} &= \int_{\Omega_e} \mathbf{B}_p^T \mathbf{s} d\Omega + \int_{\Omega_e} \nabla N^{(p)} p d\Omega \\ f_p^{\text{int},(p)} &= \int_{\Omega_e} N^{(p)} \left( \nabla \cdot \mathbf{v} - \frac{\dot{p}}{\kappa} - 3\dot{\varepsilon}^{\text{th}} \right) d\Omega \end{aligned} \right\}^{(p)} \\ \mathbf{f}_b^{\text{int}} = \int_{\Omega_e} \mathbf{B}_b^T \mathbf{s} d\Omega + \int_{\Omega_e} \nabla N^{(b)} p d\Omega \end{array} \right) \qquad (64)
$$

$$
\mathbf{f}^{\text{ext}} = \sum_{e=1}^{N_e} \left( \left\{ \begin{aligned} f_T^{\text{ext},(p)} &= \int_{\partial\Omega_e \cap \partial\Omega_q} \phi^d N^{(p)} dS + \int_{\Omega_e} N^{(p)} r d\Omega + \int_{\Omega_e} \beta \sigma : \mathbf{D}^p N^{(p)} d\Omega \\ \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 \\ 0 &= \int_{\Omega_e} \mathbf{f} N^{(b)} d\Omega \end{aligned} \right\}^{(p)} \\
\mathbf{f}_b^{\text{ext}} = \int_{\Omega_e} \mathbf{f} N^{(b)} d\Omega
$$

where the notation $\sum_{e=1}^{N_e}$ denotes the assembling operation of element quantities, $\Omega_e$ the volume of the element, and $\mathbf{B}_p$ is the matrix containing the gradients of the shape function associated to the vertex node $p$ ($1 \le p \le 4$); it is recalled that $\int_{\partial\Omega_e \cap \partial\Omega_F} \mathbf{F}_d N^{(b)} dS = \mathbf{0}$ because the bubble function vanishes on the element boundary. The heat flux vector $\mathbf{h}$ in the expression of $f_T^{\text{int},(p)}$ is related to the temperature gradient through Fourier's law:

$$
\mathbf{h} = -\lambda \nabla T. \tag{65}
$$

4.2.2. Linearization

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

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

In the same way as for the fluid-type behaviours discussed in Part I, the linearization is performed with a Newton-Raphson
method. Using an implicit Euler scheme, the following system is written at iteration *k*:

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

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)}$ is a residue defined without the inertia terms:

$$
\mathbf{R} = \mathbf{f}^{\text{ext}} - \mathbf{f}^{\text{int}}. \tag{68}
$$

The stiffness matrix $\mathbf{K}^{(k)}$ is defined as:

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

and is assembled from contributions of submatrices involving the temperature, velocity and pressure degrees of freedom of vertex nodes and the 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}). \tag{70}
$$

The system (67) is symbolically of the form:

$$
\begin{bmatrix} \mathbf{C} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix}^{(k)} \begin{Bmatrix} \dot{\mathbf{T}} \\ \dot{\mathbf{v}} \\ \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{K}_{pp} & \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 (71)
$$
where the residue is of the form:

$$ \begin{Bmatrix} \mathbf{R}_T \\ \mathbf{R}_v \\ \mathbf{R}_p \\ \mathbf{R}_b \end{Bmatrix}^{(k)} = \sum_{e=1}^{N_e} \left\{ \begin{Bmatrix} f_T^{\text{ext.}(p)} - f_T^{\text{int.}(p)} \\ \mathbf{f}_v^{\text{ext.}(p)} - \mathbf{f}_v^{\text{int.}(p)} \\ -f_p^{\text{int.}(p)} \\ \mathbf{f}_b^{\text{ext.}} - \mathbf{f}_b^{\text{int}} \end{Bmatrix} \right\}^{(k)} \quad (72) $$

Comments:

* The submatrix $K_{pp}$ appears in (71) in the case of the solid mechanics formulation, in contrast to the case of the fluid mechanics formulation presented in Part I. This submatrix results from the elastic part of the mechanical behaviour implying the weak form (61) of the internal constraint.

* The equations pertaining to the pressure degrees of freedom explicitly involve the time step through both the contribution of the elastic part of the mechanical behaviour and thermal dilation effects:

$$ (R_p^{(p)})^{(k)} = - \int_{\Omega_e} N^{(p)} \left( \nabla \cdot \mathbf{v} - \frac{1}{\kappa} \frac{\Delta p}{\Delta t} - \frac{\text{tr}\Delta\boldsymbol{\varepsilon}^{\text{th}}}{\Delta t} \right) d\Omega. \quad (73) $$

### 4.2.3. Elimination of the bubble degrees of freedom and method of solution

The degrees of freedom of the bubble node are eliminated in solid mechanics like in fluid mechanics to save CPU time. The mechanical inertia terms being neglected, the degrees of freedom of the bubble node can here be expressed explicitly as functions of those associated to vertex nodes, as done by Bellet [19]. The elimination is performed sequentially for all elements, and leads to the following system:

$$ \mathbf{C}^{(k)} \dot{\mathbf{T}}^{(k+1)} + \bar{\mathbf{K}}^{(k)} \delta \bar{\mathbf{q}}^{(k)} = \bar{\mathbf{R}}^{(k)} \quad (74) $$

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

$$ \bar{\mathbf{q}}^T = \{\mathbf{T} \cdot \mathbf{v} \cdot \mathbf{p}\}. \quad (75) $$

Afterwards, the system (74) is discretized in time leading to the following linear system at each iteration:

$$ \mathbf{A}^{(k)} \delta \bar{\mathbf{q}}^{(k)} = \mathbf{B}^{(k)} \quad (76) $$

where the matrix $\mathbf{A}^{(k)}$ and the global residue $\mathbf{B}^{(k)}$ are defined as:

$$ \mathbf{A}^{(k)} = \left( \frac{\bar{\mathbf{M}}^{(k)}}{\Delta t} + \bar{\mathbf{K}}^{(k)} \right) \quad (77) $$

$$ \mathbf{B}^{(k)} = \bar{\mathbf{R}}^{(k)} - \bar{\mathbf{M}}^{(k)} \dot{\bar{\mathbf{q}}}^{(k)}. $$

## 5. Comparison of analytical and numerical solutions

A comparison between analytical and numerical solutions is presented in this section for the elastic-plastic and thermo-elastic-plastic problems considered in Sections 2 and 3 respectively. The solid domain is discretized with 9600 P1 + /P1 (3D) finite elements and 3360 nodes; a view from above of the mesh used is shown in Fig. 5. The inner and outer radii of the mesh are respectively set to $a = 0.5$ m and $b = 1$ m.

### 5.1. Elastic-plastic problem

We are first interested in the purely mechanical problem of Section 2, therefore all thermal degrees of freedom are prescribed to zero. Thus, the numerical problem has 13 440 degrees of freedom. The displacements are set to zero at the nodes of the inner cylinder, and the displacements of the nodes of the outer cylinder are driven so that an overall rotation is prescribed. The viscometer is made of an elastic-plastic material with the following parameters:

* shear modulus: $G = 80 000$ MPa.

* yield stress in simple tension: $\sigma_0 = 200$ MPa.

* strain hardening modulus: $h = 10 000$ MPa.

An increasing loading is considered, starting from a natural initial state up to a prescribed rotation angle of 9° on the outer cylinder. The comparison between the analytical and numerical solutions is performed for rotation angles of 1°, 3°, 4°, 5°, 6° and 9°, and is shown in Fig. 6. The displacement (actually, the curved arc length swept $l$), the shear stress $\tau$ and the cumulated
<figure><img src="image_6.png" /><figcaption>Fig. 5. Mesh of the viscometer.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Comparisons of the displacement, shear stress and cumulated plastic strain.</figcaption></figure>

plastic strain $\varepsilon_{eq}$ are plotted on a radial line of the viscometer, for the six loadings defined. Discrete points correspond to numerical results and continuous lines to analytical ones. One can observe a good agreement between analytical and numerical solutions.

The evolution of the plastic crown radius can be deduced from the evolution of the cumulated plastic strain $\varepsilon_{eq}$ in Fig. 6. The geometry and the loadings have been chosen in such a way that the structure reaches the fully elastic-plastic regime for a moderate rotation of the outer cylinder. Though the solution has already reached the large strain regime (0.25 of cumulated plastic strain for the last load), the good agreement between analytical and numerical solutions confirms that the terms associated to the objective derivative of the stress disregarded in the theoretical solution are indeed negligible, at least up to the maximum loading considered here.

## 5.2. Thermo-elasto-plastic problem

In this section, a comparison of the thermal parts of the analytical and numerical thermomechanical solutions of Section 3 is presented. Since the analytical expression (54) of the temperature is applicable only when the viscometer is entirely elastic-plastic, the thermomechanical calculation is performed starting from the solution pre-computed with a purely
<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison of the temperature field.</figcaption></figure>

mechanical behaviour up to an outer rotation angle of 3°, ensuring that the viscometer is completely elastic-plastic. Since the current time is involved in the thermal part of the solution, the rotation angle is prescribed on the outer cylinder in such a way that the plastic crown radius varies exponentially in time, using Eq. (24) connecting the curved arc length swept by the outer cylinder to the plastic crown radius.

For the thermomechanical calculation, an initial condition on the temperature respecting the analytical solution (54) is prescribed on each node of the mesh. The arbitrary parameter $\mu$ arising in the analytical solution (54) is set to $50 \text{ m}^{-1}$. The numerical problem now consists of 16800 degrees of freedom as five unknowns are involved at each node of the mesh. Homogeneous Dirichlet boundary conditions are prescribed on the temperature degrees of freedom at the nodes of the inner and outer cylinders according to Eq. (52). The thermal conductivity $\lambda$, the density $\rho$ and the heat capacity $C$ are fixed to $\lambda = 600 \text{ W m}^{-1} \text{ K}^{-1}$, $\rho = 2000 \text{ kg m}^{-3}$ and $C = 150 \text{ J kg}^{-1} \text{ K}^{-1}$, independently of the temperature. Since there is no influence of the thermal part of the solution on its mechanical part, the latter part remains unchanged, see Fig. 6.

Comparisons between analytical and numerical temperatures are performed at times $t = 0.1, 0.1825, 0.2507, 0.309$ and $0.4455 \text{ s}$, corresponding respectively to rotations of the outer cylinder of $3^\circ, 4^\circ, 5^\circ, 6^\circ$ and $9^\circ$, and are plotted in Fig. 7. One can observe a very good agreement between these two solutions even for the last instant for which a very sharp peak of temperature close to the inner cylinder is apparent. This reflects the fact that the Couette viscometer thermomechanical problem for nonlinear solid-type behaviours represents a difficult numerical test involving large strains, heat generation through plastic dissipation and strong temperature gradients. The comparison just shown establishes the very good behaviour and the accuracy of the strongly coupled $P1 + /P1$ finite element in solid mechanics in a convincing manner.

## 6. Conclusion

In the present Part II, extensions of the classical Couette viscometer steady-state solution in fluid mechanics have been developed for solid-type nonlinear behaviours, first in the case of a purely quasi-static mechanical problem, and second in the case of a thermomechanical problem accounting for the generation of heat through plastic dissipation within the structure, in both small and large strains. These solutions permit to assess a new solid $P1 + /P1$ finite element developed in the framework of large strains with a temperature/velocity/pressure formulation.

In a first step, the quasi-static purely mechanical problem was addressed. The case of an elastic-plastic material obeying the von Mises criterion with isotropic hardening was considered. A solution was first found in the framework of small strains in Section 2.1. The displacement field was found by integrating a function defined from the elastic-plastic constitutive law. It was also shown that the unloading of the Couette viscometer leads to a vanishing residual stress field. The extension to the large strain framework was performed using an additive decomposition of the Eulerian strain rate, and assuming that additional terms associated to the objective derivative of the stress were negligible. The solution developed in small strain was found to be unchanged with the sole replacement of the tangential displacement by the curved arc length swept.

In a second step, the Couette viscometer problem was addressed in thermo-elasto-plasticity in Section 3, accounting for the generation of heat through plastic dissipation within the structure. The aim of the solution was to test the thermal part of the thermomechanical numerical solution computed with the $P1 + /P1$ finite element presented in Section 4. With the assumptions of material parameters independent of temperature and negligible thermal dilation effects, a solution for the thermal part of the coupled problem applicable during the fully elastic-plastic regime of the loading was found using modified Bessel's functions.
Section 4 was devoted to the description of a temperature/velocity/pressure solid P1 + /P1 finite element in the framework of large strains, associated to an implicit (backward) Euler algorithm in time. This element is dedicated to the solutions of strongly coupled thermomechanical problems arising for example in welding applications.

The assessment of this finite element was subsequently performed in Section 5. Both the mechanical and thermal parts of the numerical solution were compared to the analytical reference solutions developed. These comparisons showed a good agreement between both solutions, evidencing the good behaviour and the accuracy of the fully-coupled P1 + /P1 solid finite element developed.

## Appendix. Extension of the solution of the Couette viscometer problem in elasto-plasticity to large strains

The extension discussed is done using the Eulerian setting, involving velocities instead of displacements. The definition of the strain given by (7)₃ is replaced by that of the Eulerian strain rate:

$$2D_{r\theta} = \dot{\gamma} = \frac{dv_{\theta}}{dr} - \frac{v_{\theta}}{r} = \frac{d\dot{l}}{dr} - \frac{\dot{l}}{r} \qquad (A.1)$$

where $v_θ$ denotes the tangential component of the velocity, and $l$ the curved arc length swept by a particle. Integration yields the expression of the “slip”:

$$\gamma = \frac{dl}{dr} - \frac{l}{r} \qquad (A.2)$$

where account has been taken of the fact that $r$ remains constant in time for a given particle. Following the same approach as for small strains, the curved arc length $l$ is expressed as a function of the shear stress $\tau$ through the constitutive law:

$$\gamma = \frac{dl}{dr} - \frac{l}{r} = g(\tau(r)) \qquad (A.3)$$

where $g$ is recalled to be a function defined from the inverse function of the constitutive law. Provided that the additional terms associated to the objective derivative in the hypoelastic law are neglected, the function $g$ is still defined by Eqs. (20) and (21) in the hypoelastic and hypoelastic-plastic zones respectively. Thus the curved arc length swept by a point located at the distance $r$ from the origin is given by integration:

$$l(r) = r \int_{a}^{r} \frac{g(r')}{r'} dr'. \qquad (A.4)$$

Thus, the solution developed for small strains remains unchanged provided that the tangential displacement $u_θ$ is replaced by the curved arc length swept $l$.

## References

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

[2] 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).

[3] SYSWELD®, 2007. ESI Group. User's Manual

[4] R. Hill, The Mathematical Theory of Plasticity, Clarendon Press, Oxford, 1967.

[5] J. Salençon, De l'élastoplasticité au Calcul à la Rupture, École Polytechnique, 2002 (in French).

[6] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, second ed., Clarendon Press, Oxford, 1953.

[7] M. Abramowitz, I.A. Stegun, Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables, Courier Dover Publications, 1965.

[8] Y.G. Duan, Y. Vincent, F. Boitout, J.B. Leblond, J.M. Bergheau, Prediction of welding residual distortions of large structures using a local/global approach, J. Mech. Sci. Technol. 21 (2007) 1700-1706.

[9] S. Souloumiac, F. Boitout, J.M. Bergheau, A new local-global approach for the modelling of welded steel component distortions, Math. Model. Weld Phenom. 6 (2002) 573-590.

[10] S.A. Tsirkas, P. Papanikos, K. Pericleous, N. Strusevich, F. Boitout, J.M. Bergheau, Evaluation of distortions in laser welded shipbuilding parts using a local-global finite element approach, Sci. Technol. Weld. Joining 8 (2003) 79-88.

[11] J.M. Bergheau, Y. Vincent, J.B. Leblond, J.F. Jullien, Viscoplastic behaviour of steels during welding, Sci. Technol. Weld. Joining 9 (2004) 323-330.

[12] J.B. Leblond, D. Pont, J. Devaux, D. Bru, J.M. Bergheau, Metallurgical and mechanical consequences of phase transformations in numerical simulations of welding processes, in: Lennart Karlsson (Ed.), Modeling in Welding. Hot Powder Forming and Casting, Vol. 42, ASM International, ISBN: 0-87170-616-4, 1997, pp. 61-89. (Chapter 4).

[13] L.E. Lindgren, Finite element modeling and simulation of welding. Part I: increased complexity, J. Thermal Stresses 24 (2001) 141-192.

[14] L.E. Lindgren, Finite element modeling and simulation of welding. Part II: improved material modeling, J. Thermal Stresses 24 (2001) 195-231.

[15] L.E. Lindgren, Finite element modeling and simulation of welding. Part III: efficiency and integration, J. Thermal Stresses 24 (2001) 305-334.

[16] L.E. Lindgren, Numerical modelling of welding, Comput. Methods Appl. Mech. Eng. 195 (2006) 6710-6736.

[17] É. 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, 16-20 May 2005, 2005, pp. 261-266.

[18] É. 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. Comput. Mech. 16 (2007) 865-887 (in French).

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