# A staggered projection scheme for viscoelastic flows

Omar Mokhtari, Yohan Davit, Jean Claude Latche, Michel Quintard

► To cite this version:

Omar Mokhtari, Yohan Davit, Jean Claude Latche, Michel Quintard. A staggered projection scheme for viscoelastic flows. ESAIM: Mathematical Modelling and Numerical Analysis, 2023, 57 (3), pp.1747-1793. 10.1051/m2an/2023020. irsn-04388582

HAL Id: irsn-04388582

https://irsn.hal.science/irsn-04388582

Submitted on 11 Jan 2024

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
# A staggered projection scheme for viscoelastic flows

Omar Mokhtari, Yohan Davit, Jean Claude Latche, Michel Quintard

► To cite this version:

Omar Mokhtari, Yohan Davit, Jean Claude Latche, Michel Quintard. A staggered projection scheme for viscoelastic flows. ESAIM: Mathematical Modelling and Numerical Analysis, 2023, 57 (3), pp.1747-1793. 10.1051/m2an/2023020. irsn-04388582v2

HAL Id: irsn-04388582

https://asnr.hal.science/irsn-04388582v2

Submitted on 11 Jan 2024

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
# A STAGGERED PROJECTION SCHEME FOR VISCOELASTIC FLOWS

OMAR MOKHTARI<fnref n="1" />,<fnref n="2" />, YOHAN DAVIT<fnref n="3" />, JEAN-CLAUDE LATCHÉ<fnref n="4" />,* AND MICHEL QUINTARD<fnref n="3" />

**Abstract.** We develop a numerical scheme for the flow of viscoelastic fluids, including the Oldroyd-B and FENE-CR constitutive models. The space discretization is staggered, using either the Marker-And-Cell (MAC) scheme for structured nonuniform grids, or the Rannacher and Turek (RT) nonconforming low-order finite element approximation for general quadrangular or hexahedral meshes. The time discretization uses a fractional-step algorithm where the solution of the Navier-Stokes equations is first obtained by a projection method and then the transport-reaction equation for the conformation tensor is solved by a finite volume scheme. In order to obtain consistency, the space discretization of the divergence of the elastic part of the stress tensor in the momentum balance equation is derived using a weak form of the MAC scheme. For stability and accuracy purposes, the solution of the transport-reaction equation for the conformation tensor is split into pure convection steps, with a change of variable to the log-conformation tensor, and a reaction step, which consists in solving one ODE per cell via an Euler scheme with local sub-cycling. Numerical computations for the flow in the lid-driven cavity at Weissenberg numbers above one and the flow around a confined cylinder confirm the efficiency of the scheme.

**Mathematics Subject Classification.** 35Q31, 65N12, 65M12, 76M10, 76M12.

Received October 24, 2021. Accepted February 22, 2023.

## 1. INTRODUCTION

This paper is dedicated to presenting a new pressure-correction approach, based on a staggered arrangement of the unknowns, for viscoelastic incompressible flows of polymeric liquids. Such flows are usually modelled by solving for $(p, u, c)$ in the following coupled system of equations

$$ \rho(\partial_t u + \xi u \cdot \nabla u) = -\nabla p + \operatorname{div} \tau_s + \operatorname{div} \tau_p, \quad \text{in } (0,T) \times \Omega, \quad (1a) $$

$$ \text{with } \tau_s = \eta_s(\nabla u + (\nabla u)^t) \quad \text{and} \quad \tau_p = \frac{\eta_p}{\lambda} f(c)(c - I_d), $$

$$ \operatorname{div} u = 0, \quad \text{in } (0,T) \times \Omega, \quad (1b) $$

$$ \partial_t c + u \cdot \nabla c = (\nabla u) c + c (\nabla u)^t - \frac{1}{\lambda} g(c)(c - I_d), \quad \text{in } (0,T) \times \Omega, \quad (1c) $$

*Keywords and phrases.* Viscoelastic flows, finite volume, staggered discretization, projection scheme.

<footnote n="1">¹ Institut National Polytechnique de Toulouse (INPT), Toulouse, France.</footnote>

<footnote n="2">² TotalEnergies E&P, CSTJF, Pau, France.</footnote>

<footnote n="3">³ Institut de Mécanique des Fluides de Toulouse (IMFT), CNRS & Université de Toulouse, Toulouse, France.</footnote>

<footnote n="4">⁴ Institut de Radioprotection et de Sûreté Nucléaire (IRSN), Fontenay-aux-Roses, France.</footnote>

<footnote>*Corresponding author: jean-claude.latche@irsn.fr</footnote>
# A STAGGERED PROJECTION SCHEME FOR VISCOELASTIC FLOWS

OMAR MOKHTARI<fnref n="1" />,<fnref n="2" />, YOHAN DAVIT<fnref n="3" />, JEAN-CLAUDE LATCHÉ<fnref n="4" />,* AND MICHEL QUINTARD<fnref n="3" />

**Abstract.** We develop a numerical scheme for the flow of viscoelastic fluids, including the Oldroyd-B and FENE-CR constitutive models. The space discretization is staggered, using either the Marker-And-Cell (MAC) scheme for structured nonuniform grids, or the Rannacher and Turek (RT) nonconforming low-order finite element approximation for general quadrangular or hexahedral meshes. The time discretization uses a fractional-step algorithm where the solution of the Navier-Stokes equations is first obtained by a projection method and then the transport-reaction equation for the conformation tensor is solved by a finite volume scheme. In order to obtain consistency, the space discretization of the divergence of the elastic part of the stress tensor in the momentum balance equation is derived using a weak form of the MAC scheme. For stability and accuracy purposes, the solution of the transport-reaction equation for the conformation tensor is split into pure convection steps, with a change of variable to the log-conformation tensor, and a reaction step, which consists in solving one ODE per cell via an Euler scheme with local sub-cycling. Numerical computations for the flow in the lid-driven cavity at Weissenberg numbers above one and the flow around a confined cylinder confirm the efficiency of the scheme.

**Mathematics Subject Classification.** 35Q31, 65N12, 65M12, 76M10, 76M12.

Received October 24, 2021. Accepted February 22, 2023.

## 1. INTRODUCTION

This paper is dedicated to presenting a new pressure-correction approach, based on a staggered arrangement of the unknowns, for viscoelastic incompressible flows of polymeric liquids. Such flows are usually modelled by solving for $(p, u, c)$ in the following coupled system of equations

$$ \rho(\partial_t u + \xi u \cdot \nabla u) = -\nabla p + \text{div } \tau_s + \text{div } \tau_p, \quad \text{in } (0,T) \times \Omega, \quad (1a) $$

$$ \text{with } \tau_s = \eta_s(\nabla u + (\nabla u)^t) \quad \text{and} \quad \tau_p = \frac{\eta_p}{\lambda} f(c)(c - I_d), $$

$$ \text{div } u = 0, \quad \text{in } (0,T) \times \Omega, \quad (1b) $$

$$ \partial_t c + u \cdot \nabla c = (\nabla u) c + c (\nabla u)^t - \frac{1}{\lambda} g(c)(c - I_d), \quad \text{in } (0,T) \times \Omega, \quad (1c) $$

*Keywords and phrases.* Viscoelastic flows, finite volume, staggered discretization, projection scheme.

<footnote n="1">¹ Institut National Polytechnique de Toulouse (INPT), Toulouse, France.</footnote>

<footnote n="2">² TotalEnergies E&P, CSTJF, Pau, France.</footnote>

<footnote n="3">³ Institut de Mécanique des Fluides de Toulouse (IMFT), CNRS & Université de Toulouse, Toulouse, France.</footnote>

<footnote n="4">⁴ Institut de Radioprotection et de Sûreté Nucléaire (IRSN), Fontenay-aux-Roses, France.</footnote>

<footnote>*Corresponding author: jean-claude.latche@irsn.fr</footnote>
where $\Omega$ is an open bounded connected subset of $\mathbb{R}^d$ with $d \in \{2, 3\}$, $(0, T)$ a finite time interval with $T > 0$, $\mathbf{u}$ the velocity field of the fluid, $p$ the pressure, and $\mathbf{c}$ the symmetric positive-definite polymer conformation tensor. The operator $\mathbf{u} \cdot \nabla$ is defined by

$$ (\mathbf{u} \cdot \nabla \mathbf{a})_i = \mathbf{u} \cdot \nabla a_i = \operatorname{div}(a_i \mathbf{u}) \quad \text{or} \quad (\mathbf{u} \cdot \nabla \mathbf{a})_{i,j} = \mathbf{u} \cdot \nabla a_{i,j} = \operatorname{div}(a_{i,j} \mathbf{u}), \quad 1 \le i, j \le d, $$

for $\mathbf{a}$ a generic vector or tensor field, respectively. The constant coefficients $\eta_s$, $\eta_p$, $\rho$, and $\lambda$ are respectively the solvent viscosity, the polymer viscosity, the fluid density and the polymer retardation time, and $I_d$ stands for the identity matrix of $\mathbb{R}^{d \times d}$. The coefficient $\xi$ serves to artificially eliminate $\mathbf{u} \cdot \nabla \mathbf{u}$, so as to obtain the unsteady Stokes equations for $\xi = 0$ and the Navier-Stokes equations for $\xi = 1$. The tensor $\tau_s$ is the part of the stress accounting for the fluid solvent and $\tau_p$ is the part of the stress accounting for the presence of polymers. The functions $f(\mathbf{c})$ and $g(\mathbf{c})$ vary depending on the viscoelastic constitutive equations (see [10, 45] for a review) with, for example, the Oldroyd-B model given by $f(\mathbf{c}) = g(\mathbf{c}) = I_d$ and the FENE-CR model by $f(\mathbf{c}) = g(\mathbf{c}) = \frac{b}{b-\operatorname{tr}(\mathbf{c})} I_d$, with $b$ a real number greater than the space dimension $d$. The dimensionless parameters that characterize the flow are the Reynolds number, $Re = \rho UL/(\eta_s + \eta_p)$, and the Weissenberg number, $Wi = \lambda U/L$, where $U$ and $L$ are a characteristic velocity and length scale. This system must be complemented by initial conditions for the velocity and the conformation tensor and by suitable boundary conditions.

A variety of numerical schemes have been proposed for the resolution of this system (see [2] for a recent review). All of these schemes face the long-standing high-Weissenberg number problem (HWNP), which makes numerical simulations extremely difficult as soon as the Weissenberg number gets close to 1. In practice, steep gradients of the stress field rapidly develop, making it difficult to correctly capture the profile of the conformation tensor and to guarantee its positivity. One of the most significant advances in pushing back the limits of the HWNP was the log-formulation proposed in [22, 23] and used in many works [18, 20, 42, 43, 46, 47]. It consists in reformulating the problem using a logarithmic transformation of the conformation tensor. The idea described by the authors is that the conformation tensor varies exponentially and it is therefore easier to approach the logarithm of the conformation tensor with polynomial approximations than the tensor itself. Further, this formulation naturally maintains the positive-definiteness of the conformation tensor. It has later given birth to a variety of similar approaches, such as the square-root transformation [4] or the kernel conformation transformation [1]. Despite these advances, the HWNP also often limits applications to moderate Weissenberg numbers because of computational costs and much progress is still needed.

One of the most important challenges of the HWNP can be understood by considering the one-dimensional transport of the scalar $c$ posed, say, on $(0, L_1) \times (0, T)$,

$$ \partial_t c + \partial_x c = a(x)c, \quad c(x, 0) = c_0(x), \quad x \in (0, L_1) \text{ and } c(0, t) = c_i, \quad t \in (0, T), \quad (2) $$

where $a(x)$ is here to mimick the velocity gradients. This problem is similar to that of [23] but with the difference that we consider here that $a(x)$ can become very large and sharply varies in space. When this is the case, obtaining accurate solutions of this toy problem and, by extension, of the HWNP becomes particularly difficult. This toy problem actually belongs to a class of hyperbolic conservation laws with stiff source terms [17, 38, 39] that are found, for example, in modelling shallow-water flows with non-flat bottoms [3, 21]. A classical method, such as a centered discretization of the source term, yields non-physical solutions [9]. The primary difficulty is to preserve the steady-state solutions in the discrete version of the operators. Schemes that verify this property are usually termed “well-balanced”, as introduced in [11, 26, 33]. One of the simplest ways to do this is to incorporate the discrete equilibrium states into the fluxes, which is the principle of the Upwinding Sources at Interfaces (“U.S.I.”) schemes [9, 11, 12, 34]. These methods are much more accurate than classical approaches but seem next to impossible to generalise to non-scalar 2D or 3D problems. Some viscoelastic flows that admit a steady-state solution (for instance, the so-called lid-driven cavity problem, at least at moderate Wi numbers) may be schematically understood by considering the steady transport problem $\partial_x c = a(x)c$ with a field $a(x)$ generating, by a change of sign, successively opposite variations of $c$. This may be understood as variations of
where $\Omega$ is an open bounded connected subset of $\mathbb{R}^d$ with $d \in \{2, 3\}$, $(0, T)$ a finite time interval with $T > 0$, $\mathbf{u}$ the velocity field of the fluid, $p$ the pressure, and $\mathbf{c}$ the symmetric positive-definite polymer conformation tensor. The operator $\mathbf{u} \cdot \nabla$ is defined by

$$ (\mathbf{u} \cdot \nabla \mathbf{a})_i = \mathbf{u} \cdot \nabla a_i = \operatorname{div}(a_i \mathbf{u}) \quad \text{or} \quad (\mathbf{u} \cdot \nabla \mathbf{a})_{i,j} = \mathbf{u} \cdot \nabla a_{i,j} = \operatorname{div}(a_{i,j} \mathbf{u}), \quad 1 \le i, j \le d, $$

for $\mathbf{a}$ a generic vector or tensor field, respectively. The constant coefficients $\eta_s$, $\eta_p$, $\rho$, and $\lambda$ are respectively the solvent viscosity, the polymer viscosity, the fluid density and the polymer retardation time, and $I_d$ stands for the identity matrix of $\mathbb{R}^{d \times d}$. The coefficient $\xi$ serves to artificially eliminate $\mathbf{u} \cdot \nabla \mathbf{u}$, so as to obtain the unsteady Stokes equations for $\xi = 0$ and the Navier-Stokes equations for $\xi = 1$. The tensor $\tau_s$ is the part of the stress accounting for the fluid solvent and $\tau_p$ is the part of the stress accounting for the presence of polymers. The functions $f(\mathbf{c})$ and $g(\mathbf{c})$ vary depending on the viscoelastic constitutive equations (see [10, 45] for a review) with, for example, the Oldroyd-B model given by $f(\mathbf{c}) = g(\mathbf{c}) = I_d$ and the FENE-CR model by $f(\mathbf{c}) = g(\mathbf{c}) = \frac{b}{b-\operatorname{tr}(\mathbf{c})} I_d$, with $b$ a real number greater than the space dimension $d$. The dimensionless parameters that characterize the flow are the Reynolds number, $Re = \rho UL/(\eta_s + \eta_p)$, and the Weissenberg number, $Wi = \lambda U/L$, where $U$ and $L$ are a characteristic velocity and length scale. This system must be complemented by initial conditions for the velocity and the conformation tensor and by suitable boundary conditions.

A variety of numerical schemes have been proposed for the resolution of this system (see [2] for a recent review). All of these schemes face the long-standing high-Weissenberg number problem (HWNP), which makes numerical simulations extremely difficult as soon as the Weissenberg number gets close to 1. In practice, steep gradients of the stress field rapidly develop, making it difficult to correctly capture the profile of the conformation tensor and to guarantee its positivity. One of the most significant advances in pushing back the limits of the HWNP was the log-formulation proposed in [22, 23] and used in many works [18, 20, 42, 43, 46, 47]. It consists in reformulating the problem using a logarithmic transformation of the conformation tensor. The idea described by the authors is that the conformation tensor varies exponentially and it is therefore easier to approach the logarithm of the conformation tensor with polynomial approximations than the tensor itself. Further, this formulation naturally maintains the positive-definiteness of the conformation tensor. It has later given birth to a variety of similar approaches, such as the square-root transformation [4] or the kernel conformation transformation [1]. Despite these advances, the HWNP also often limits applications to moderate Weissenberg numbers because of computational costs and much progress is still needed.

One of the most important challenges of the HWNP can be understood by considering the one-dimensional transport of the scalar $c$ posed, say, on $(0, L_1) \times (0, T)$,

$$ \partial_t c + \partial_x c = a(x)c, \quad c(x, 0) = c_0(x), \quad x \in (0, L_1) \text{ and } c(0, t) = c_i, \quad t \in (0, T), \quad (2) $$

where $a(x)$ is here to mimick the velocity gradients. This problem is similar to that of [23] but with the difference that we consider here that $a(x)$ can become very large and sharply varies in space. When this is the case, obtaining accurate solutions of this toy problem and, by extension, of the HWNP becomes particularly difficult. This toy problem actually belongs to a class of hyperbolic conservation laws with stiff source terms [17, 38, 39] that are found, for example, in modelling shallow-water flows with non-flat bottoms [3, 21]. A classical method, such as a centered discretization of the source term, yields non-physical solutions [9]. The primary difficulty is to preserve the steady-state solutions in the discrete version of the operators. Schemes that verify this property are usually termed “well-balanced”, as introduced in [11, 26, 33]. One of the simplest ways to do this is to incorporate the discrete equilibrium states into the fluxes, which is the principle of the Upwinding Sources at Interfaces (“U.S.I.”) schemes [9, 11, 12, 34]. These methods are much more accurate than classical approaches but seem next to impossible to generalise to non-scalar 2D or 3D problems. Some viscoelastic flows that admit a steady-state solution (for instance, the so-called lid-driven cavity problem, at least at moderate Wi numbers) may be schematically understood by considering the steady transport problem $\partial_x c = a(x)c$ with a field $a(x)$ generating, by a change of sign, successively opposite variations of $c$. This may be understood as variations of
the tensor $c$ along closed streamlines. In Appendix D, we solve this steady-state problem with large values of $a(x)$ to mimic the HWNP. We show that a logarithmic transformation of the transport equation drastically improves the accuracy of the scheme, even though it does not preserve the steady-state and is therefore not well-balanced.

It has also been suggested in [37, 40] that a numerical scheme for viscoelastic flows should preserve, at the discrete level, a *free-energy* estimate derived at the continuum-level. For the system of equations (1), the time derivative of the Helmholtz *free-energy* is equal to the opposite of the dissipation terms plus energy sources from volume forcing terms or boundary conditions. For an initial-value problem with no sources and homogeneous Dirichlet boundary conditions, the Helmholtz *free-energy* must therefore decrease in time (see [7, 8, 31]). An ideal numerical scheme should preserve this behaviour at the discrete level and, in so doing, ensure the long-time stability of solutions. Following this idea, Barrett and Boyaval [5, 6] and Boyaval *et al.* [14] developed numerical schemes that verify the *free-energy* estimate and proved the global-in-time existence of solutions to those discretizations, together with convergence to weak solutions of the continuous problem, under time-step restriction or with an additional dissipative term in the conformation equation. The authors also suggest in [13] (see also [41]) that this approach is not sufficient to guarantee accuracy.

Here, we aim for a numerical scheme that

(1) preserves the positive definiteness of the tensor $c$,

(2) satisfies a discrete equivalent of the *free-energy* estimate,

(3) minimizes CPU costs and is well suited to high-performance computations,

(4) shows a good accuracy, in particular to capture steady-state solutions,

(5) can be readily used for a variety of viscoelastic constitutive laws.

We present in Sections 2 and 3 an algorithm that satisfies Point (1) and, for the Oldroyd-B model, Point (2). We use a piecewise constant space approximation for the conformation tensor, which *in fine* (and especially for the final algorithm, see next paragraph) offers great advantages from an algorithmic point of view. Equation (1c) is then solved using a finite-volume scheme. To ensure the so-called *inf-sup* stability while keeping a consistent order of accuracy (and so, of CPU cost), we combine this discretization with a staggered-in-space approximation of the velocity/pressure pair, with either the Marker-And-Cell (MAC) scheme [28, 29], or a nonconforming low-order finite element, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedric cells (extension to simplicial cells with the Crouzeix-Raviart element [19] is straightforward). To keep the algorithm efficient and robust, we introduce a fractional-step algorithm where the solution of the Navier-Stokes equations (1a) and (1b) is computed by a standard projection method. Satisfying the discrete *free-energy* estimate has however a price: the transport equation for $c$ must remain coupled with the correction step, with the end-of-step conformation tensor and velocity used in both this transport equation and the velocity correction equation. This yields a set of two coupled nonlinear equations for the pressure and $c$, after using the usual trick that consists in combining the velocity correction and the divergence constraint into an elliptic problem for the pressure only. Finally, satisfying the *free-energy* estimate also requires a discrete weak formulation of (1a) and, in particular, of the diffusion term and of $\operatorname{div} \boldsymbol{\tau}_p$. This is in fact natural for finite element discretizations, but not so for the MAC scheme for which this discretization is thus derived in Section 3.3. We also observe and prove that this weak formulation enjoys a built-in Lax-Wendroff weak consistency property [24].

Unfortunately, respecting Points (3) and (4) requires that we make a compromise regarding the time discretization, as detailed in Section 4. We proceed by decoupling the system of equations for the projection step, using an explicit-in-time approximation of the conformation tensor in the velocity correction equation, which greatly reduces the cost of each iteration but breaks the proof of the *free-energy* estimate. Consequently, Point (2) is not strictly satisfied anymore: the space semi-discretization of the momentum balance equation is consistent with the *free-energy* estimate, but not the time discretization. For accuracy, we further introduce a log transformation of the (now decoupled) transport equation for the conformation tensor. Unlike the standard approach of [23], however, this transformation is only partial. In the spirit of [43], we proceed by operator splitting in equation (1c), with a change of variable from $c$ to $\log(c)$ in pure convection steps and a separate
the tensor $c$ along closed streamlines. In Appendix D, we solve this steady-state problem with large values of $a(x)$ to mimic the HWNP. We show that a logarithmic transformation of the transport equation drastically improves the accuracy of the scheme, even though it does not preserve the steady-state and is therefore not well-balanced.

It has also been suggested in [37, 40] that a numerical scheme for viscoelastic flows should preserve, at the discrete level, a *free-energy* estimate derived at the continuum-level. For the system of equations (1), the time derivative of the Helmholtz *free-energy* is equal to the opposite of the dissipation terms plus energy sources from volume forcing terms or boundary conditions. For an initial-value problem with no sources and homogeneous Dirichlet boundary conditions, the Helmholtz *free-energy* must therefore decrease in time (see [7, 8, 31]). An ideal numerical scheme should preserve this behaviour at the discrete level and, in so doing, ensure the long-time stability of solutions. Following this idea, Barrett and Boyaval [5, 6] and Boyaval *et al.* [14] developed numerical schemes that verify the *free-energy* estimate and proved the global-in-time existence of solutions to those discretizations, together with convergence to weak solutions of the continuous problem, under time-step restriction or with an additional dissipative term in the conformation equation. The authors also suggest in [13] (see also [41]) that this approach is not sufficient to guarantee accuracy.

Here, we aim for a numerical scheme that

(1) preserves the positive definiteness of the tensor $c$,

(2) satisfies a discrete equivalent of the *free-energy* estimate,

(3) minimizes CPU costs and is well suited to high-performance computations,

(4) shows a good accuracy, in particular to capture steady-state solutions,

(5) can be readily used for a variety of viscoelastic constitutive laws.

We present in Sections 2 and 3 an algorithm that satisfies Point (1) and, for the Oldroyd-B model, Point (2). We use a piecewise constant space approximation for the conformation tensor, which *in fine* (and especially for the final algorithm, see next paragraph) offers great advantages from an algorithmic point of view. Equation (1c) is then solved using a finite-volume scheme. To ensure the so-called *inf-sup* stability while keeping a consistent order of accuracy (and so, of CPU cost), we combine this discretization with a staggered-in-space approximation of the velocity/pressure pair, with either the Marker-And-Cell (MAC) scheme [28, 29], or a nonconforming low-order finite element, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedric cells (extension to simplicial cells with the Crouzeix-Raviart element [19] is straightforward). To keep the algorithm efficient and robust, we introduce a fractional-step algorithm where the solution of the Navier-Stokes equations (1a) and (1b) is computed by a standard projection method. Satisfying the discrete *free-energy* estimate has however a price: the transport equation for $c$ must remain coupled with the correction step, with the end-of-step conformation tensor and velocity used in both this transport equation and the velocity correction equation. This yields a set of two coupled nonlinear equations for the pressure and $c$, after using the usual trick that consists in combining the velocity correction and the divergence constraint into an elliptic problem for the pressure only. Finally, satisfying the *free-energy* estimate also requires a discrete weak formulation of (1a) and, in particular, of the diffusion term and of $\operatorname{div} \tau_p$. This is in fact natural for finite element discretizations, but not so for the MAC scheme for which this discretization is thus derived in Section 3.3. We also observe and prove that this weak formulation enjoys a built-in Lax-Wendroff weak consistency property [24].

Unfortunately, respecting Points (3) and (4) requires that we make a compromise regarding the time discretization, as detailed in Section 4. We proceed by decoupling the system of equations for the projection step, using an explicit-in-time approximation of the conformation tensor in the velocity correction equation, which greatly reduces the cost of each iteration but breaks the proof of the *free-energy* estimate. Consequently, Point (2) is not strictly satisfied anymore: the space semi-discretization of the momentum balance equation is consistent with the *free-energy* estimate, but not the time discretization. For accuracy, we further introduce a log transformation of the (now decoupled) transport equation for the conformation tensor. Unlike the standard approach of [23], however, this transformation is only partial. In the spirit of [43], we proceed by operator splitting in equation (1c), with a change of variable from $c$ to $\log(c)$ in pure convection steps and a separate
reaction step without any transformation. This latter step consists in solving one ODE per mesh cell thanks to the piecewise constant discretization of $c$ and, in contrast with [43], these ODEs are solved directly for $c$, and not $\log(c)$, so as to avoid any artificial introduction of nonlinearities, together with a cumbersome factorization of the velocity gradient tensor. The viscoelastic constitutive law is only involved in this cell-by-cell procedure (up to an explicit computation of $\tau_p$ for the velocity prediction step), so that it can be easily modified to capture a variety of viscoelastic behaviors (Point (5)). To further gain in efficiency, we use a local time step for the solution of the ODE on each mesh element; this allows us to adjust the time step for each cell to a value small enough to preserve the positive definiteness of $c$ (so Point (1) is still satisfied, see Appendix A.2 for the proof) while preventing a blow-up of the CPU cost. Indeed, in practice, the velocity gradient is large only in a small number of cells where the ODE becomes stiff and subcycling is needed. This greatly limits the computational cost, although one must be careful in the domain decomposition for parallel computations to keep the computing charge of the processors equilibrated. Summing up, our approach is thus based upon a low-order numerical scheme that is very efficient and reasonably accurate, the most important cause of inaccuracy being tackled by the log transformation in the transport step and a time-step adaptation to cope with the stiffness of the reaction terms. The constitutive laws for the viscoelastic fluids is essentially isolated in the resolution of cell-by-cell ODEs and can be easily modified. The resulting solver is available and open-source, as a specific module of the CALIF³S platform [16] developed at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN).

In the last part of the paper, Section 5, we study the efficiency and accuracy of the scheme using two benchmark 2D tests, the lid-driven cavity (up to a Weissenberg of 5) and the cylinder confined in a channel. Using massively parallel computations, we obtain solutions for the cylinder problem that are, to our knowledge, more accurate than already published results. In addition, we take benefit of the versatility of our algorithm with respect to the viscoelastic constitutive law to perform a systematic comparison of the results obtained with the Oldroyd-B and the FENE-CR models.

## 2. MESHES AND UNKNOWNS

Let $\mathcal{M}$ be a decomposition of the domain $\Omega$, supposed to be regular in the usual sense of the finite element literature. The cells may be:

- quadrilaterals ($d=2$) or hexahedra ($d=3$) for a general polygonal or polyhedral domain $\Omega$,

- rectangles ($d=2$) or rectangular parallelepipeds ($d=3$) for a domain $\Omega$ with boundaries that are hyperplanes normal to a coordinate axis. In this case, the mesh is defined by a union of Cartesian non-uniform grids.

The sets $\mathcal{E}$ and $\mathcal{E}(K)$ are, respectively, the set of all $(d-1)$-faces $\sigma$ of the mesh and of the element $K \in \mathcal{M}$. The set of faces included in the boundary of $\Omega$ is $\mathcal{E}_{\text{ext}}$ and the set of internal faces is $\mathcal{E}_{\text{int}} = \mathcal{E} \setminus \mathcal{E}_{\text{ext}}$. A face $\sigma \in \mathcal{E}_{\text{int}}$ separating the cells $K$ and $L$ is denoted by $\sigma = K|L$. The vector $n_{K,\sigma}$ is the unit outward normal vector to a face $\sigma$ of $K$. For $1 \le i \le d$, we denote by $\mathcal{E}^{(i)}$, $\mathcal{E}_{\text{int}}^{(i)}$ and $\mathcal{E}_{\text{ext}}^{(i)}$ the subset of the faces of $\mathcal{E}$, $\mathcal{E}_{\text{int}}$ and $\mathcal{E}_{\text{ext}}$, respectively, that are perpendicular to the $i$th unit vector of the canonical basis of $\mathbb{R}^d$ (these definitions being useful only in the case of Cartesian meshes). For $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}$, $|K|$ is the measure of $K$ and $|\sigma|$ the $(d-1)$-measure of the face $\sigma$.

The space discretization is staggered, using either the Marker-And-Cell (MAC) scheme [28, 29] or non-conforming low-order finite element approximations, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedral meshes. The degrees of freedom for the pressure and the conformation tensor (i.e. the discrete pressure and conformation tensor unknowns) are associated with the cells of the mesh $\mathcal{M}$ and are denoted by $\{p_K, c_K, K \in \mathcal{M}\}$. For the velocity, we have the following cases:

- For the Rannacher-Turek discretization, the discrete velocity unknowns are located at the center of the faces of the mesh and we choose the version of the element where they represent the average of the velocity through a face. For the sake of simplicity, we suppose in this presentation that the velocity is prescribed to zero over the whole boundary of $\Omega$ and these Dirichlet boundary conditions are taken into account by setting
reaction step without any transformation. This latter step consists in solving one ODE per mesh cell thanks to the piecewise constant discretization of $c$ and, in contrast with [43], these ODEs are solved directly for $c$, and not $\log(c)$, so as to avoid any artificial introduction of nonlinearities, together with a cumbersome factorization of the velocity gradient tensor. The viscoelastic constitutive law is only involved in this cell-by-cell procedure (up to an explicit computation of $\tau_p$ for the velocity prediction step), so that it can be easily modified to capture a variety of viscoelastic behaviors (Point (5)). To further gain in efficiency, we use a local time step for the solution of the ODE on each mesh element; this allows us to adjust the time step for each cell to a value small enough to preserve the positive definiteness of $c$ (so Point (1) is still satisfied, see Appendix A.2 for the proof) while preventing a blow-up of the CPU cost. Indeed, in practice, the velocity gradient is large only in a small number of cells where the ODE becomes stiff and subcycling is needed. This greatly limits the computational cost, although one must be careful in the domain decomposition for parallel computations to keep the computing charge of the processors equilibrated. Summing up, our approach is thus based upon a low-order numerical scheme that is very efficient and reasonably accurate, the most important cause of inaccuracy being tackled by the log transformation in the transport step and a time-step adaptation to cope with the stiffness of the reaction terms. The constitutive laws for the viscoelastic fluids is essentially isolated in the resolution of cell-by-cell ODEs and can be easily modified. The resulting solver is available and open-source, as a specific module of the CALIF³S platform [16] developed at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN).

In the last part of the paper, Section 5, we study the efficiency and accuracy of the scheme using two benchmark 2D tests, the lid-driven cavity (up to a Weissenberg of 5) and the cylinder confined in a channel. Using massively parallel computations, we obtain solutions for the cylinder problem that are, to our knowledge, more accurate than already published results. In addition, we take benefit of the versatility of our algorithm with respect to the viscoelastic constitutive law to perform a systematic comparison of the results obtained with the Oldroyd-B and the FENE-CR models.

## 2. MESHES AND UNKNOWNS

Let $\mathcal{M}$ be a decomposition of the domain $\Omega$, supposed to be regular in the usual sense of the finite element literature. The cells may be:

- quadrilaterals ($d=2$) or hexahedra ($d=3$) for a general polygonal or polyhedral domain $\Omega$,

- rectangles ($d=2$) or rectangular parallelepipeds ($d=3$) for a domain $\Omega$ with boundaries that are hyperplanes normal to a coordinate axis. In this case, the mesh is defined by a union of Cartesian non-uniform grids.

The sets $\mathcal{E}$ and $\mathcal{E}(K)$ are, respectively, the set of all $(d-1)$-faces $\sigma$ of the mesh and of the element $K \in \mathcal{M}$. The set of faces included in the boundary of $\Omega$ is $\mathcal{E}_{\text{ext}}$ and the set of internal faces is $\mathcal{E}_{\text{int}} = \mathcal{E} \setminus \mathcal{E}_{\text{ext}}$. A face $\sigma \in \mathcal{E}_{\text{int}}$ separating the cells $K$ and $L$ is denoted by $\sigma = K|L$. The vector $n_{K,\sigma}$ is the unit outward normal vector to a face $\sigma$ of $K$. For $1 \le i \le d$, we denote by $\mathcal{E}^{(i)}$, $\mathcal{E}_{\text{int}}^{(i)}$ and $\mathcal{E}_{\text{ext}}^{(i)}$ the subset of the faces of $\mathcal{E}$, $\mathcal{E}_{\text{int}}$ and $\mathcal{E}_{\text{ext}}$, respectively, that are perpendicular to the $i$th unit vector of the canonical basis of $\mathbb{R}^d$ (these definitions being useful only in the case of Cartesian meshes). For $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}$, $|K|$ is the measure of $K$ and $|\sigma|$ the $(d-1)$-measure of the face $\sigma$.

The space discretization is staggered, using either the Marker-And-Cell (MAC) scheme [28, 29] or non-conforming low-order finite element approximations, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedral meshes. The degrees of freedom for the pressure and the conformation tensor (i.e. the discrete pressure and conformation tensor unknowns) are associated with the cells of the mesh $\mathcal{M}$ and are denoted by $\{p_K, c_K, K \in \mathcal{M}\}$. For the velocity, we have the following cases:

- For the Rannacher-Turek discretization, the discrete velocity unknowns are located at the center of the faces of the mesh and we choose the version of the element where they represent the average of the velocity through a face. For the sake of simplicity, we suppose in this presentation that the velocity is prescribed to zero over the whole boundary of $\Omega$ and these Dirichlet boundary conditions are taken into account by setting
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Primal and dual cells for the Rannacher–Turek elements.</figcaption></figure>

the velocity unknowns associated with an external face to zero. The set of discrete velocity unknowns reads {$u_{\sigma,i}, \sigma \in \mathcal{E}_{\text{int}}, 1 \le i \le d$}.

- For the MAC discretization, the degrees of freedom for the ith component of the velocity are located at the centre of the faces $\sigma \in \mathcal{E}^{(i)}$, so, with the same boundary conditions, the whole set of discrete velocity unknowns reads {$u_{\sigma,i}, \sigma \in \mathcal{E}_{\text{int}}^{(i)}, 1 \le i \le d$}.

Hence, there are d unknowns per face of the primal mesh in the case of the RT scheme, corresponding to the d components of the velocity, and there is only one unknown per face of the primal mesh in the case of the MAC scheme, corresponding to the normal component of the velocity.

We now introduce a dual mesh for the finite volume approximation of the time derivative and convection terms in the momentum balance equation.

- For the Rannacher–Turek discretization, the dual mesh is the same for all the velocity components. When $K \in \mathcal{M}$ is a rectangle or a cuboid, we define $D_{K,\sigma}$, for $\sigma \in \mathcal{E}(K)$, as the cone with basis $\sigma$ and with vertex the mass center of $K$ (see Fig. 1). We thus obtain a partition of $K$ in $m$ sub-volumes, where $m$ is the number of faces of the cell ($m = 4$ in two space dimensions and $m = 6$ in three space dimensions), each sub-volume having the same measure $|D_{K,\sigma}| = |K|/m$. We extend this definition to general quadrangles and hexahedra by supposing that we have built a partition still of equal-volume sub-cells and with the same connectivities. The volume $D_{K,\sigma}$ is referred to as the half-diamond cell associated with $K$ and $\sigma$. Note that the shape of $D_{K,\sigma}$ does not need to be specified (and, for highly distorted meshes, if $|\sigma|$ is too small, $D_{K,\sigma}$ cannot be a cone). For $\sigma \in \mathcal{E}_{\text{int}}$, $\sigma = K|L$, we now define the diamond cell $D_{\sigma}$ associated with $\sigma$ by $D_{\sigma} = D_{K,\sigma} \cup D_{L,\sigma}$.

- For the MAC discretization, the dual mesh depends on the component of the velocity. For each component, the MAC dual mesh only differs from the RT dual mesh by the choice of the half-diamond cell, which, for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}(K)$, is now the rectangle or rectangular parallelepiped of basis $\sigma$ and of measure $|D_{K,\sigma}| = |K|/2$ (see Fig. 2).

We denote by $|D_{\sigma}|$ the measure of the dual cell $D_{\sigma}$ and by $\varepsilon = D_{\sigma}|D_{\sigma'}$ the face separating two diamond cells $D_{\sigma}$ and $D_{\sigma'}$. The set of the (dual) faces of $D_{\sigma}$ is denoted by $\bar{\mathcal{E}}(D_{\sigma})$.

Finally, in order to be able to write a unique expression of the discrete equations for both MAC and RT scheme, we introduce the set of faces $\mathcal{E}_S^{(i)}$ associated with the degrees of freedom of the ith component of the velocity (S stands for “scheme”):

$$ \mathcal{E}_S^{(i)} = \begin{cases} \mathcal{E}_{\text{int}}^{(i)} & \text{for the MAC scheme,} \\ \mathcal{E}_{\text{int}} & \text{for the RT scheme.} \end{cases} $$
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Primal and dual cells for the Rannacher–Turek elements.</figcaption></figure>

the velocity unknowns associated with an external face to zero. The set of discrete velocity unknowns reads {$u_{\sigma,i}, \sigma \in \mathcal{E}_{\text{int}}, 1 \le i \le d$}.

- For the MAC discretization, the degrees of freedom for the ith component of the velocity are located at the centre of the faces $\sigma \in \mathcal{E}^{(i)}$, so, with the same boundary conditions, the whole set of discrete velocity unknowns reads {$u_{\sigma,i}, \sigma \in \mathcal{E}_{\text{int}}^{(i)}, 1 \le i \le d$}.

Hence, there are d unknowns per face of the primal mesh in the case of the RT scheme, corresponding to the d components of the velocity, and there is only one unknown per face of the primal mesh in the case of the MAC scheme, corresponding to the normal component of the velocity.

We now introduce a dual mesh for the finite volume approximation of the time derivative and convection terms in the momentum balance equation.

- For the Rannacher–Turek discretization, the dual mesh is the same for all the velocity components. When $K \in \mathcal{M}$ is a rectangle or a cuboid, we define $D_{K,\sigma}$, for $\sigma \in \mathcal{E}(K)$, as the cone with basis $\sigma$ and with vertex the mass center of $K$ (see Fig. 1). We thus obtain a partition of $K$ in $m$ sub-volumes, where $m$ is the number of faces of the cell ($m = 4$ in two space dimensions and $m = 6$ in three space dimensions), each sub-volume having the same measure $|D_{K,\sigma}| = |K|/m$. We extend this definition to general quadrangles and hexahedra by supposing that we have built a partition still of equal-volume sub-cells and with the same connectivities. The volume $D_{K,\sigma}$ is referred to as the half-diamond cell associated with $K$ and $\sigma$. Note that the shape of $D_{K,\sigma}$ does not need to be specified (and, for highly distorted meshes, if $|\sigma|$ is too small, $D_{K,\sigma}$ cannot be a cone). For $\sigma \in \mathcal{E}_{\text{int}}$, $\sigma = K|L$, we now define the diamond cell $D_{\sigma}$ associated with $\sigma$ by $D_{\sigma} = D_{K,\sigma} \cup D_{L,\sigma}$.

- For the MAC discretization, the dual mesh depends on the component of the velocity. For each component, the MAC dual mesh only differs from the RT dual mesh by the choice of the half-diamond cell, which, for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}(K)$, is now the rectangle or rectangular parallelepiped of basis $\sigma$ and of measure $|D_{K,\sigma}| = |K|/2$ (see Fig. 2).

We denote by $|D_{\sigma}|$ the measure of the dual cell $D_{\sigma}$ and by $\varepsilon = D_{\sigma}|D_{\sigma'}$ the face separating two diamond cells $D_{\sigma}$ and $D_{\sigma'}$. The set of the (dual) faces of $D_{\sigma}$ is denoted by $\bar{\mathcal{E}}(D_{\sigma})$.

Finally, in order to be able to write a unique expression of the discrete equations for both MAC and RT scheme, we introduce the set of faces $\mathcal{E}_S^{(i)}$ associated with the degrees of freedom of the ith component of the velocity (S stands for “scheme”):

$$ \mathcal{E}_S^{(i)} = \begin{cases} \mathcal{E}_{\text{int}}^{(i)} & \text{for the MAC scheme,} \\ \mathcal{E}_{\text{int}} & \text{for the RT scheme.} \end{cases} $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Primal and dual cells for the MAC discretization in two space dimensions. (a) Mesh associated to the first component of the velocity. (b) Mesh associated to the second component of the velocity.</figcaption></figure>

### 3. A FREE-ENERGY PRESERVING PRESSURE CORRECTION SCHEME FOR OLDROYD-B FLUIDS

Let us first recall the formal derivation of the energy estimate given in [31]. Here, we restrict the exposition to the Oldroyd-B model, and the starting point is system (1) where we set $f(c) = g(c) = I_d$. We consider that no-slip conditions are prescribed on the whole boundary. Consequently, there is no source of energy in the system, either through volume forces or boundary conditions, so that we expect a decrease of the free-energy due to dissipation processes. The first step to show this consists in taking the inner product of the momentum balance equation (1a) with $u$, which yields

$$ \partial_t E_{\text{kin}} + \xi \mathbf{u} \cdot \nabla E_{\text{kin}} = \left( -\nabla p + \operatorname{div} \tau_s(\mathbf{u}) + \operatorname{div} \left( \frac{\eta_p}{\lambda} \mathbf{c} \right) \right) \cdot \mathbf{u}, \quad (3) $$

where $E_{\text{kin}}$ is the kinetic energy defined as

$$ E_{\text{kin}} = \frac{1}{2} \rho |\mathbf{u}|^2. $$

Let us also consider the trace of equation (1c) multiplied by $\eta_p/(2\lambda)$. In doing so, we suppose that any tensor $c$ solution of equation (1c) is symmetric definite positive. We have the relationships $\operatorname{tr}(\nabla u c) = \nabla u : c$ and $\operatorname{tr}(c(\nabla u)^t) = \nabla u : c$, thanks to the fact that $\operatorname{tr}(c) = \operatorname{tr}(c^t)$ and $\operatorname{tr}(cd) = c : d^t$ for any symmetric matrix $c$ of $\mathbb{R}^{d \times d}$ and any matrix $d$ of $\mathbb{R}^{d \times d}$. We obtain

$$ \frac{\eta_p}{2\lambda} (\partial_t \operatorname{tr}(c) + \mathbf{u} \cdot \nabla \operatorname{tr}(c)) = \frac{\eta_p}{\lambda} c : \nabla \mathbf{u} - \frac{\eta_p}{2\lambda^2} \operatorname{tr}(c - I_d). \quad (4) $$

Since $c$ is supposed to be invertible, we may also multiply equation (1c) by $c^{-1}$ to get

$$ c^{-1} \partial_t c + c^{-1} \mathbf{u} \cdot \nabla c = c^{-1} (\nabla \mathbf{u}) c + (\nabla \mathbf{u})^t - \frac{1}{\lambda} (I_d - c^{-1}). \quad (5) $$
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Primal and dual cells for the MAC discretization in two space dimensions. (a) Mesh associated to the first component of the velocity. (b) Mesh associated to the second component of the velocity.</figcaption></figure>

### 3. A FREE-ENERGY PRESERVING PRESSURE CORRECTION SCHEME FOR OLDROYD-B FLUIDS

Let us first recall the formal derivation of the energy estimate given in [31]. Here, we restrict the exposition to the Oldroyd-B model, and the starting point is system (1) where we set $f(c) = g(c) = I_d$. We consider that no-slip conditions are prescribed on the whole boundary. Consequently, there is no source of energy in the system, either through volume forces or boundary conditions, so that we expect a decrease of the free-energy due to dissipation processes. The first step to show this consists in taking the inner product of the momentum balance equation (1a) with $u$, which yields

$$ \partial_t E_{\text{kin}} + \xi \mathbf{u} \cdot \nabla E_{\text{kin}} = \left( -\nabla p + \operatorname{div} \tau_s(\mathbf{u}) + \operatorname{div} \left( \frac{\eta_p}{\lambda} \mathbf{c} \right) \right) \cdot \mathbf{u}, \quad (3) $$

where $E_{\text{kin}}$ is the kinetic energy defined as

$$ E_{\text{kin}} = \frac{1}{2} \rho |\mathbf{u}|^2. $$

Let us also consider the trace of equation (1c) multiplied by $\eta_p/(2\lambda)$. In doing so, we suppose that any tensor $c$ solution of equation (1c) is symmetric definite positive. We have the relationships $\operatorname{tr}(\nabla u c) = \nabla u : c$ and $\operatorname{tr}(c(\nabla u)^t) = \nabla u : c$, thanks to the fact that $\operatorname{tr}(c) = \operatorname{tr}(c^t)$ and $\operatorname{tr}(cd) = c : d^t$ for any symmetric matrix $c$ of $\mathbb{R}^{d \times d}$ and any matrix $d$ of $\mathbb{R}^{d \times d}$. We obtain

$$ \frac{\eta_p}{2\lambda} \left( \partial_t \operatorname{tr}(c) + \mathbf{u} \cdot \nabla \operatorname{tr}(c) \right) = \frac{\eta_p}{\lambda} c : \nabla \mathbf{u} - \frac{\eta_p}{2\lambda^2} \operatorname{tr}(c - I_d). \quad (4) $$

Since $c$ is supposed to be invertible, we may also multiply equation (1c) by $c^{-1}$ to get

$$ c^{-1} \partial_t c + c^{-1} \mathbf{u} \cdot \nabla c = c^{-1} (\nabla \mathbf{u}) c + (\nabla \mathbf{u})^t - \frac{1}{\lambda} (I_d - c^{-1}). \quad (5) $$
We are now going to consider the trace of this equation. For the first term on the left-hand side, we know (see [14], Lem. 1.2) that tr($c^{-1}\partial_t c$) = $\partial_t$tr(ln(c)). For the second term on the left-hand side, consider that

$$u \cdot \nabla c = \sum_{k=1}^{d} u_k \partial_k c,$$

with $\partial_k c$ the tensor obtained by deriving each entry of c with respect to the kth space coordinate. By the same identity as for the first term, we therefore have tr($c^{-1}u \cdot \nabla c$) = $u \cdot \nabla$tr(ln(c)). For the first term on the right-hand side, we remark that tr($c^{-1} \nabla u c$) = tr(cc^{-1} \nabla u) = tr(\nabla u) = div(u) = 0$ since trcd) = tr(d c) for any matrix c and d of $\mathbb{R}^{d \times d}$. For the second term on the right-hand side, we also have tr($\nabla u)^t$) = 0. Therefore, taking the trace of equation (5) and multiplying by $\frac{\eta_p}{2\lambda}$, we get

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(\ln(c)) + u \cdot \nabla \mathrm{tr}(\ln(c))) = -\frac{\eta_p}{2\lambda^2} \mathrm{tr}(I_d - c^{-1}), \quad (6)$$

which can also be written as

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(\ln(c) + I_d) + u \cdot \nabla \mathrm{tr}(\ln(c) + I_d)) = -\frac{\eta_p}{2\lambda^2} \mathrm{tr}(I_d - c^{-1}), \quad (7)$$

considering that

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(I_d) + u \cdot \nabla \mathrm{tr}(I_d)) = 0. \quad (8)$$

We can now build a transport equation for the free-energy by adding equation (3) to equation (4) and subtracting equation (7) to obtain

$$\partial_t E + u \cdot \nabla (\xi E_{\text{kin}} + E_p) = (-\nabla p + \text{div } \tau_s(u) + \text{div} \left( \frac{\eta_p}{\lambda} c \right)) \cdot u + \frac{\eta_p}{\lambda} c : \nabla u - \frac{\eta_p}{2\lambda^2} \text{tr}(c + c^{-1} - 2I_d),$$

where E, the free-energy, and $E_p$ are defined by

$$E = E_{\text{kin}} + E_p, \quad E_p = \frac{\eta_p}{2\lambda} \mathrm{tr}(c - \ln(c) - I_d).$$

Finally, integrating over $\Omega$ yields

$$\frac{d}{dt} \int_{\Omega} E \, dx + \int_{\Omega} \tau_s(u) : \nabla u \, dx + \frac{\eta_p}{2\lambda^2} \int_{\Omega} \mathrm{tr}(c + c^{-1} - 2I_d) \, dx = 0, \quad (9)$$

by integration by parts and using the no-slip boundary condition. For any symmetric positive definite tensor c, both c - ln(c) - I_d and c + c^{-1} - 2I_d are symmetric positive definite ([14], Lem. 1.1). Thus, tr(c + c^{-1} - 2I_d) > 0, tr(c - ln(c) - I_d) > 0 and the previous relation shows that the free-energy E decreases with time, since the viscous dissipation $\tau_s(u) : \nabla u$ is non-negative. Equation (9) thus provides a stability estimate for the solutions of (1).

### 3.1. The time semi-discrete scheme

Our aim is now to derive a time semi-discretization of system (1), still for the Oldroyd-B model, satisfying a semi-discrete analogue of equation (9). Let a uniform partition $0 = t_0 < t_1 < \dots < t_N = T$ of $(0, T)$ be given, with a constant time step $\delta t$. We consider the following semi-discretization in time of System (1), which belongs to the class of pressure correction schemes.

**Prediction step** – Solve for $\tilde{u}^{n+1}$:

$$\frac{\rho}{\delta t} (\tilde{u}^{n+1} - u^n) + \xi \rho u^n \cdot \nabla \tilde{u}^{n+1} - \mathrm{div}(\tau_s(\tilde{u}^{n+1})) + \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, c^n = 0. \quad (10a)$$
We are now going to consider the trace of this equation. For the first term on the left-hand side, we know (see [14], Lem. 1.2) that tr($c^{-1}\partial_t c$) = $\partial_t$tr(ln(c)). For the second term on the left-hand side, consider that

$$u \cdot \nabla c = \sum_{k=1}^{d} u_k \partial_k c,$$

with $\partial_k c$ the tensor obtained by deriving each entry of c with respect to the kth space coordinate. By the same identity as for the first term, we therefore have tr($c^{-1}u \cdot \nabla c$) = $u \cdot \nabla$tr(ln(c)). For the first term on the right-hand side, we remark that tr($c^{-1}\nabla u c$) = tr(cc^{-1}\nabla u) = tr(\nabla u) = div(u) = 0$ since trcd) = tr(d c) for any matrix c and d of $\mathbb{R}^{d \times d}$. For the second term on the right-hand side, we also have tr($\nabla u)^t$) = 0. Therefore, taking the trace of equation (5) and multiplying by $\frac{\eta_p}{2\lambda}$, we get

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(\ln(c)) + u \cdot \nabla \mathrm{tr}(\ln(c))) = -\frac{\eta_p}{2\lambda^2} \mathrm{tr}(I_d - c^{-1}), \quad (6)$$

which can also be written as

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(\ln(c) + I_d) + u \cdot \nabla \mathrm{tr}(\ln(c) + I_d)) = -\frac{\eta_p}{2\lambda^2} \mathrm{tr}(I_d - c^{-1}), \quad (7)$$

considering that

$$\frac{\eta_p}{2\lambda} (\partial_t \mathrm{tr}(I_d) + u \cdot \nabla \mathrm{tr}(I_d)) = 0. \quad (8)$$

We can now build a transport equation for the free-energy by adding equation (3) to equation (4) and subtracting equation (7) to obtain

$$\partial_t E + u \cdot \nabla (\xi E_{\text{kin}} + E_p) = (-\nabla p + \text{div } \tau_s(u) + \text{div} \left( \frac{\eta_p}{\lambda} c \right)) \cdot u + \frac{\eta_p}{\lambda} c : \nabla u - \frac{\eta_p}{2\lambda^2} \text{tr}(c + c^{-1} - 2I_d),$$

where E, the free-energy, and $E_p$ are defined by

$$E = E_{\text{kin}} + E_p, \quad E_p = \frac{\eta_p}{2\lambda} \mathrm{tr}(c - \ln(c) - I_d).$$

Finally, integrating over $\Omega$ yields

$$\frac{d}{dt} \int_{\Omega} E \, dx + \int_{\Omega} \tau_s(u) : \nabla u \, dx + \frac{\eta_p}{2\lambda^2} \int_{\Omega} \mathrm{tr}(c + c^{-1} - 2I_d) \, dx = 0, \quad (9)$$

by integration by parts and using the no-slip boundary condition. For any symmetric positive definite tensor c, both c - ln(c) - I_d and c + c^{-1} - 2I_d are symmetric positive definite ([14], Lem. 1.1). Thus, tr(c + c^{-1} - 2I_d) > 0, tr(c - ln(c) - I_d) > 0 and the previous relation shows that the free-energy E decreases with time, since the viscous dissipation $\tau_s(u) : \nabla u$ is non-negative. Equation (9) thus provides a stability estimate for the solutions of (1).

## 3.1. The time semi-discrete scheme

Our aim is now to derive a time semi-discretization of system (1), still for the Oldroyd-B model, satisfying a semi-discrete analogue of equation (9). Let a uniform partition $0 = t_0 < t_1 < \dots < t_N = T$ of $(0, T)$ be given, with a constant time step $\delta t$. We consider the following semi-discretization in time of System (1), which belongs to the class of pressure correction schemes.

**Prediction step** – Solve for $\tilde{u}^{n+1}$:

$$\frac{\rho}{\delta t} (\tilde{u}^{n+1} - u^n) + \xi \rho u^n \cdot \nabla \tilde{u}^{n+1} - \mathrm{div}(\tau_s(\tilde{u}^{n+1})) + \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, c^n = 0. \quad (10a)$$
**Correction step** – Solve for $p^{n+1}$ and $\mathbf{u}^{n+1}$:

$$ \frac{\rho}{\delta t} (\mathbf{u}^{n+1} - \tilde{\mathbf{u}}^{n+1}) + \nabla (p^{n+1} - p^n) - \frac{\eta_p}{\lambda} \mathrm{div} (\mathbf{c}^{n+1} - \mathbf{c}^n) = 0, \quad (10b) $$

$$ \mathrm{div} \, \mathbf{u}^{n+1} = 0, \quad (10c) $$

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1} - \mathbf{c}^n) + \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} = \nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1} + \mathbf{c}^{n+1} (\nabla \mathbf{u}^{n+1})^t - \frac{1}{\lambda} (\mathbf{c}^{n+1} - \mathbf{I}_d). \quad (10d) $$

This time semi-discretization satisfies the stability estimate stated in the following result.

**Proposition 3.1.** Let $(\mathbf{u}^n, p^n, \mathbf{c}^n)_{0\le n\le N}$ be the solution of System (10). For $0 \le n \le N$, let us suppose that $\mathbf{c}^n$ is positive definite and let the discrete free-energy of the solution be defined by

$$ E^n(\mathbf{u}^n, p^n, \mathbf{c}^n, \delta t) = \frac{\rho}{2} |\mathbf{u}^n|^2 + \frac{\eta_p}{2\lambda} \mathrm{tr}(\mathbf{c}^n - \ln(\mathbf{c}^n) - \mathbf{I}_d) + \frac{\delta t^2}{2\rho} |\nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n|^2. \quad (11) $$

Then we have the following stability estimate

$$ \int_{\Omega} E^{n+1} \mathrm{d}\mathbf{x} + \frac{\rho}{2} \int_{\Omega} |\tilde{\mathbf{u}} - \mathbf{u}^n|^2 \mathrm{d}\mathbf{x} + \delta t \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x} \\ + \frac{\eta_p}{2\lambda^2} \int_{\Omega} \mathrm{tr}(\mathbf{c}^{n+1} - (\mathbf{c}^{n+1})^{-1} - 2\mathbf{I}_d) \mathrm{d}\mathbf{x} \le \int_{\Omega} E^n \mathrm{d}\mathbf{x}. \quad (12) $$

Consequently, the integral over the domain of the discrete free-energy is non-increasing.

*Proof.* Using the identity $2a(a-b) = a^2 + (a-b)^2 - b^2$ valid for any pair of real numbers a and b, we have

$$ \tilde{\mathbf{u}}^{n+1} \cdot \left( \frac{\rho}{\delta t} (\tilde{\mathbf{u}}^{n+1} - \mathbf{u}^n) + \xi \rho \mathbf{u}^n \cdot \nabla \tilde{\mathbf{u}}^{n+1} \right) \ge \frac{\rho}{2\delta t} \left( |\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2 \right) + \xi \frac{\rho}{2} \mathbf{u}^n \cdot \nabla |\tilde{\mathbf{u}}^{n+1}|^2. $$

Hence, using the fact that $\mathbf{u}^n$ is divergence-free to recast the transport term as a convection term, thanks to the no-slip boundary condition,

$$ \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \left( \frac{\rho}{\delta t} (\tilde{\mathbf{u}}^{n+1} - \mathbf{u}^n) + \xi \rho \mathbf{u}^n \cdot \nabla \tilde{\mathbf{u}}^{n+1} \right) \mathrm{d}\mathbf{x} \ge \int_{\Omega} \frac{\rho}{2\delta t} \left( |\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2 \right) \mathrm{d}\mathbf{x}. \quad (13) $$

Thanks once again to the boundary conditions, we have

$$ - \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \mathrm{div}(\tau_s(\tilde{\mathbf{u}}^{n+1})) \mathrm{d}\mathbf{x} = \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x}. \quad (14) $$

Therefore, taking the inner product of equation (10a) with $\tilde{\mathbf{u}}^{n+1}$ and integrating over the computational domain yields

$$ \frac{\rho}{2\delta t} \int_{\Omega} (|\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2) \mathrm{d}\mathbf{x} + \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x} + \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \left( \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n \right) \mathrm{d}\mathbf{x} \le 0. \quad (15) $$

Let us now rewrite equation (10b) as

$$ \frac{\rho}{\delta t} \mathbf{u}^{n+1} + \nabla p^{n+1} - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^{n+1} = \frac{\rho}{\delta t} \tilde{\mathbf{u}}^{n+1} + \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n. $$

Taking the square of the $L^2$-norm of both sides and multiplying by $\delta t/(2\rho)$ yields

$$ \frac{\rho}{2\delta t} \int_{\Omega} |\mathbf{u}^{n+1}|^2 \mathrm{d}\mathbf{x} - \int_{\Omega} \frac{\eta_p}{\lambda} \mathbf{u}^{n+1} \cdot \mathrm{div} \, \mathbf{c}^{n+1} \mathrm{d}\mathbf{x} + \frac{\delta t}{2\rho} \int_{\Omega} |\nabla p^{n+1} - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^{n+1}|^2 \mathrm{d}\mathbf{x} $$
**Correction step** – Solve for $p^{n+1}$ and $\mathbf{u}^{n+1}$:

$$ \frac{\rho}{\delta t} (\mathbf{u}^{n+1} - \tilde{\mathbf{u}}^{n+1}) + \nabla (p^{n+1} - p^n) - \frac{\eta_p}{\lambda} \mathrm{div} (\mathbf{c}^{n+1} - \mathbf{c}^n) = 0, \quad (10b) $$

$$ \mathrm{div} \, \mathbf{u}^{n+1} = 0, \quad (10c) $$

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1} - \mathbf{c}^n) + \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} = \nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1} + \mathbf{c}^{n+1} (\nabla \mathbf{u}^{n+1})^t - \frac{1}{\lambda} (\mathbf{c}^{n+1} - \mathbf{I}_d). \quad (10d) $$

This time semi-discretization satisfies the stability estimate stated in the following result.

**Proposition 3.1.** Let $(\mathbf{u}^n, p^n, \mathbf{c}^n)_{0\le n\le N}$ be the solution of System (10). For $0 \le n \le N$, let us suppose that $\mathbf{c}^n$ is positive definite and let the discrete free-energy of the solution be defined by

$$ E^n(\mathbf{u}^n, p^n, \mathbf{c}^n, \delta t) = \frac{\rho}{2} |\mathbf{u}^n|^2 + \frac{\eta_p}{2\lambda} \mathrm{tr}(\mathbf{c}^n - \ln(\mathbf{c}^n) - \mathbf{I}_d) + \frac{\delta t^2}{2\rho} |\nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n|^2. \quad (11) $$

Then we have the following stability estimate

$$ \int_{\Omega} E^{n+1} \mathrm{d}\mathbf{x} + \frac{\rho}{2} \int_{\Omega} |\tilde{\mathbf{u}} - \mathbf{u}^n|^2 \mathrm{d}\mathbf{x} + \delta t \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x} \\ + \frac{\eta_p}{2\lambda^2} \int_{\Omega} \mathrm{tr}(\mathbf{c}^{n+1} - (\mathbf{c}^{n+1})^{-1} - 2\mathbf{I}_d) \mathrm{d}\mathbf{x} \le \int_{\Omega} E^n \mathrm{d}\mathbf{x}. \quad (12) $$

Consequently, the integral over the domain of the discrete free-energy is non-increasing.

*Proof.* Using the identity $2a(a-b) = a^2 + (a-b)^2 - b^2$ valid for any pair of real numbers a and b, we have

$$ \tilde{\mathbf{u}}^{n+1} \cdot \left( \frac{\rho}{\delta t} (\tilde{\mathbf{u}}^{n+1} - \mathbf{u}^n) + \xi \rho \mathbf{u}^n \cdot \nabla \tilde{\mathbf{u}}^{n+1} \right) \ge \frac{\rho}{2\delta t} \left( |\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2 \right) + \xi \frac{\rho}{2} \mathbf{u}^n \cdot \nabla |\tilde{\mathbf{u}}^{n+1}|^2. $$

Hence, using the fact that $\mathbf{u}^n$ is divergence-free to recast the transport term as a convection term, thanks to the no-slip boundary condition,

$$ \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \left( \frac{\rho}{\delta t} (\tilde{\mathbf{u}}^{n+1} - \mathbf{u}^n) + \xi \rho \mathbf{u}^n \cdot \nabla \tilde{\mathbf{u}}^{n+1} \right) \mathrm{d}\mathbf{x} \ge \int_{\Omega} \frac{\rho}{2\delta t} \left( |\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2 \right) \mathrm{d}\mathbf{x}. \quad (13) $$

Thanks once again to the boundary conditions, we have

$$ - \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \mathrm{div}(\tau_s(\tilde{\mathbf{u}}^{n+1})) \mathrm{d}\mathbf{x} = \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x}. \quad (14) $$

Therefore, taking the inner product of equation (10a) with $\tilde{\mathbf{u}}^{n+1}$ and integrating over the computational domain yields

$$ \frac{\rho}{2\delta t} \int_{\Omega} (|\tilde{\mathbf{u}}^{n+1}|^2 - |\mathbf{u}^n|^2) \mathrm{d}\mathbf{x} + \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} \mathrm{d}\mathbf{x} + \int_{\Omega} \tilde{\mathbf{u}}^{n+1} \cdot \left( \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n \right) \mathrm{d}\mathbf{x} \le 0. \quad (15) $$

Let us now rewrite equation (10b) as

$$ \frac{\rho}{\delta t} \mathbf{u}^{n+1} + \nabla p^{n+1} - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^{n+1} = \frac{\rho}{\delta t} \tilde{\mathbf{u}}^{n+1} + \nabla p^n - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^n. $$

Taking the square of the $L^2$-norm of both sides and multiplying by $\delta t/(2\rho)$ yields

$$ \frac{\rho}{2\delta t} \int_{\Omega} |\mathbf{u}^{n+1}|^2 \mathrm{d}\mathbf{x} - \int_{\Omega} \frac{\eta_p}{\lambda} \mathbf{u}^{n+1} \cdot \mathrm{div} \, \mathbf{c}^{n+1} \mathrm{d}\mathbf{x} + \frac{\delta t}{2\rho} \int_{\Omega} |\nabla p^{n+1} - \frac{\eta_p}{\lambda} \mathrm{div} \, \mathbf{c}^{n+1}|^2 \mathrm{d}\mathbf{x} $$
$$ = \frac{\rho}{2\delta t} \int_{\Omega} |\tilde{u}^{n+1}|^2 d\mathbf{x} + \int_{\Omega} \tilde{u}^{n+1} \cdot \left( \nabla p^n - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^n \right) d\mathbf{x} + \frac{\delta t}{2\rho} \int_{\Omega} \left| \nabla p^n - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^n \right|^2 d\mathbf{x}. $$

Note that, to obtain this relation, we have used the following velocity divergence-pressure gradient duality (which holds thanks to the no-slip boundary conditions)

$$ \int_{\Omega} \nabla p^{n+1} \cdot \mathbf{u}^{n+1} d\mathbf{x} = - \int_{\Omega} p^{n+1} \operatorname{div}(\mathbf{u}^{n+1}) d\mathbf{x} = 0. \quad (16) $$

Adding with (15), we obtain the following discrete kinetic energy inequality

$$ \frac{\rho}{2\delta t} \int_{\Omega} (|\mathbf{u}^{n+1}|^2 - |\mathbf{u}^n|^2) d\mathbf{x} + \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} d\mathbf{x} - \int_{\Omega} \frac{\eta_p}{\lambda} \mathbf{u}^{n+1} \cdot \operatorname{div} \mathbf{c}^{n+1} d\mathbf{x} + T^{n+1} - T^n \le 0, \quad (17) $$

with, for $k = n$ and $k = n + 1$,

$$ T^k = \frac{\delta t}{2\rho} \int_{\Omega} \left| \nabla p^k - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^k \right|^2 d\mathbf{x}. $$

Let us now multiply (10d) by $(\mathbf{c}^{n+1})^{-1}$ to obtain

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1})^{-1} (\mathbf{c}^{n+1} - \mathbf{c}^n) + (\mathbf{c}^{n+1})^{-1} \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} = (\mathbf{c}^{n+1})^{-1} \nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1} + (\nabla \mathbf{u}^{n+1})^t - \frac{1}{\lambda} (\mathbf{I}_d - (\mathbf{c}^{n+1})^{-1}). \quad (18) $$

For any symmetric definite positive matrices $c$ and $d$, we have $\operatorname{tr}(c^{-1}(c-d)) \le \operatorname{tr}\ln(c) - \operatorname{tr}\ln(d)$ (see [14], Lem. 2.1), which yields

$$ \frac{1}{\delta t} \left( \operatorname{tr} \ln(\mathbf{c}^{n+1}) - \operatorname{tr} \ln(\mathbf{c}^n) \right) + \mathbf{u}^n \cdot \nabla \operatorname{tr} \ln(\mathbf{c}^{n+1}) \ge \operatorname{tr} \left( \frac{1}{\delta t} (\mathbf{c}^{n+1})^{-1} (\mathbf{c}^{n+1} - \mathbf{c}^n) + (\mathbf{c}^{n+1})^{-1} \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} \right). \quad (19) $$

Taking the trace of relation (18) and integrating over $\Omega$, by arguments already used in the continuous setting, we thus get

$$ \int_{\Omega} \operatorname{tr} \ln(\mathbf{c}^{n+1}) - \operatorname{tr} \ln(\mathbf{c}^n) d\mathbf{x} \ge - \int_{\Omega} \frac{1}{\lambda} \operatorname{tr}(\mathbf{I}_d - (\mathbf{c}^{n+1})^{-1}) d\mathbf{x}. $$

Taking the trace of equation (10d), adding the discrete kinetic energy balance (17) and subtracting the latter relation as in the continuous setting yields the desired estimate, thanks to the duality relation

$$ \int_{\Omega} \operatorname{tr}(\nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1}) d\mathbf{x} = \int_{\Omega} \nabla \mathbf{u}^{n+1} : \mathbf{c}^{n+1} d\mathbf{x} = - \int_{\Omega} \mathbf{u}^{n+1} \cdot \operatorname{div} \mathbf{c}^{n+1} d\mathbf{x}. \quad (20) $$

## 3.2. Towards a fully discrete scheme

A discretization of the scheme (10) is obtained by combining a finite volume and a finite element approach for the approximation of the momentum balance equation and using a finite volume technique for the transport of the conformation tensor. The general form of the scheme reads:

**Prediction step** – Solve for $\tilde{\mathbf{u}}^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, $$

$$ \frac{\rho}{\delta t} (\tilde{u}_{\sigma,i}^{n+1} - u_{\sigma,i}^n) + \xi \rho \operatorname{div}_{\sigma} (\tilde{u}_i^{n+1} \mathbf{u}^n) - \operatorname{div}_{\sigma,i} \tau_s (\tilde{u}^{n+1}) + \nabla_{\sigma,i} p^n - \operatorname{div}_{\sigma,i} \tau_p^n = 0. \quad (21a) $$
$$ = \frac{\rho}{2\delta t} \int_{\Omega} |\tilde{u}^{n+1}|^2 d\mathbf{x} + \int_{\Omega} \tilde{u}^{n+1} \cdot \left( \nabla p^n - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^n \right) d\mathbf{x} + \frac{\delta t}{2\rho} \int_{\Omega} \left| \nabla p^n - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^n \right|^2 d\mathbf{x}. $$

Note that, to obtain this relation, we have used the following velocity divergence-pressure gradient duality (which holds thanks to the no-slip boundary conditions)

$$ \int_{\Omega} \nabla p^{n+1} \cdot \mathbf{u}^{n+1} d\mathbf{x} = - \int_{\Omega} p^{n+1} \operatorname{div}(\mathbf{u}^{n+1}) d\mathbf{x} = 0. \quad (16) $$

Adding with (15), we obtain the following discrete kinetic energy inequality

$$ \frac{\rho}{2\delta t} \int_{\Omega} (|\mathbf{u}^{n+1}|^2 - |\mathbf{u}^n|^2) d\mathbf{x} + \int_{\Omega} \tau_s(\tilde{\mathbf{u}}^{n+1}) : \nabla \tilde{\mathbf{u}}^{n+1} d\mathbf{x} - \int_{\Omega} \frac{\eta_p}{\lambda} \mathbf{u}^{n+1} \cdot \operatorname{div} \mathbf{c}^{n+1} d\mathbf{x} + T^{n+1} - T^n \le 0, \quad (17) $$

with, for $k=n$ and $k=n+1$,

$$ T^k = \frac{\delta t}{2\rho} \int_{\Omega} \left| \nabla p^k - \frac{\eta_p}{\lambda} \operatorname{div} \mathbf{c}^k \right|^2 d\mathbf{x}. $$

Let us now multiply (10d) by $(\mathbf{c}^{n+1})^{-1}$ to obtain

$$ \frac{1}{\delta t} (\mathbf{c}^{n+1})^{-1} (\mathbf{c}^{n+1} - \mathbf{c}^n) + (\mathbf{c}^{n+1})^{-1} \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} = (\mathbf{c}^{n+1})^{-1} \nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1} + (\nabla \mathbf{u}^{n+1})^t - \frac{1}{\lambda} (\mathbf{I}_d - (\mathbf{c}^{n+1})^{-1}). \quad (18) $$

For any symmetric definite positive matrices $c$ and $d$, we have $\operatorname{tr}(c^{-1}(c-d)) \le \operatorname{tr}\ln(c) - \operatorname{tr}\ln(d)$ (see [14], Lem. 2.1), which yields

$$ \frac{1}{\delta t} \left( \operatorname{tr} \ln(\mathbf{c}^{n+1}) - \operatorname{tr} \ln(\mathbf{c}^n) \right) + \mathbf{u}^n \cdot \nabla \operatorname{tr} \ln(\mathbf{c}^{n+1}) \ge \operatorname{tr} \left( \frac{1}{\delta t} (\mathbf{c}^{n+1})^{-1} (\mathbf{c}^{n+1} - \mathbf{c}^n) + (\mathbf{c}^{n+1})^{-1} \mathbf{u}^n \cdot \nabla \mathbf{c}^{n+1} \right). \quad (19) $$

Taking the trace of relation (18) and integrating over $\Omega$, by arguments already used in the continuous setting, we thus get

$$ \int_{\Omega} \operatorname{tr} \ln(\mathbf{c}^{n+1}) - \operatorname{tr} \ln(\mathbf{c}^n) d\mathbf{x} \ge - \int_{\Omega} \frac{1}{\lambda} \operatorname{tr}(\mathbf{I}_d - (\mathbf{c}^{n+1})^{-1}) d\mathbf{x}. $$

Taking the trace of equation (10d), adding the discrete kinetic energy balance (17) and subtracting the latter relation as in the continuous setting yields the desired estimate, thanks to the duality relation

$$ \int_{\Omega} \operatorname{tr}(\nabla \mathbf{u}^{n+1} \mathbf{c}^{n+1}) d\mathbf{x} = \int_{\Omega} \nabla \mathbf{u}^{n+1} : \mathbf{c}^{n+1} d\mathbf{x} = - \int_{\Omega} \mathbf{u}^{n+1} \cdot \operatorname{div} \mathbf{c}^{n+1} d\mathbf{x}. \quad (20) $$

## 3.2. Towards a fully discrete scheme

A discretization of the scheme (10) is obtained by combining a finite volume and a finite element approach for the approximation of the momentum balance equation and using a finite volume technique for the transport of the conformation tensor. The general form of the scheme reads:

**Prediction step** – Solve for $\tilde{\mathbf{u}}^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, $$

$$ \frac{\rho}{\delta t} (\tilde{u}_{\sigma,i}^{n+1} - u_{\sigma,i}^n) + \xi \rho \operatorname{div}_{\sigma} (\tilde{u}_i^{n+1} \mathbf{u}^n) - \operatorname{div}_{\sigma,i} \tau_s (\tilde{u}^{n+1}) + \nabla_{\sigma,i} p^n - \operatorname{div}_{\sigma,i} \tau_p^n = 0. \quad (21a) $$
**Correction step** – Solve for $p^{n+1}$, $u^{n+1}$ and $c^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \qquad \frac{\rho}{\delta t} (u_{\sigma,i}^{n+1} - \tilde{u}_{\sigma,i}^{n+1}) + \nabla_{\sigma,i} (p^{n+1} - p^n) - \text{div}_{\sigma,i} \tau_p^{n+1} + \text{div}_{\sigma,i} \tau_p^n = 0, \qquad (21b) $$

$$ \forall K \in \mathcal{M}, \quad \text{div}_K(u^{n+1}) = 0, \qquad (21c) $$

$$ \forall K \in \mathcal{M}, \qquad \frac{1}{\delta t} (c_K^{n+1} - c_K^n) + \text{div}_K(u^n c^{n+1}) = \nabla_K u^{n+1} c_K^{n+1} + (c_K^{n+1})^t (\nabla_K u^{n+1})^t - \frac{1}{\lambda} (c_K^{n+1} - I_d). \quad (21d) $$

In the prediction step and in the velocity correction step, for $k = n$ and $k = n + 1$, the tensor $\tau_p^k$ is discretized on the primal cells and is computed as a function of the conformation tensor as

$$ (\tau_p^k)_K = \frac{\eta_p}{\lambda} (c_K^k - I_d), \quad \text{for } K \in \mathcal{M}. $$

The divergence operator applied to this tensor is defined in Section 3.3.

To build a fully discrete scheme that satisfies the same stability estimate as its time-discrete counterpart, we need to consider all the properties that have been used in the proof of Proposition 3.1 and satisfy them at the discrete level by a suitable definition of the operators in the scheme. First, we require that the transport equation for the tensor $c$ preserves its symmetry and its positive-definiteness. The former property is easy to satisfy by solving a “symmetrized version” of the right-hand side (i.e. a discretization of $\nabla u c + c^t (\nabla u)^t$ instead of $\nabla u c + c (\nabla u)^t$), as is done in equation (21d). Once the equation is symmetrized, if $c^n$ is symmetric, for $1 \le i, j \le d$, the discrete equation for $c_{i,j}^{n+1}$ is strictly the same as the equation for $c_{j,i}^{n+1}$. Conversely, a symmetric tensor solution to (21d) is also a solution of the original equation and thus switching to (21d) preserves the consistency of the scheme. Moreover, it is proved in Appendix A that a first-order backward upwind scheme for (21d) preserves the positive definiteness of $c$, under a condition for the time-step depending on the $L^\infty$-norm of $\nabla u^{n+1}$. Finally, for this equation and this convection scheme, we need a discrete analogue of (19) which is proved in Appendix C.

Let us now turn to the momentum balance equation. A finite-volume discretization of the convection term for the RT discretization may be found, in the more general context of variable density flows, in Section 4.2 of [35]. The discrete kinetic energy balance (i.e. the discrete equivalent of (13)) is established in Lemma 7.2 of [35]. For the MAC scheme, a discretization of the convection term is proposed in Section 2, equation (29) of [24] and its $L^2$-stability (i.e. as before, the discrete equivalent of (13)) is shown in Lemma 3.6 of [24]. The convection operator takes the following general form

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \quad \rho \text{div}_\sigma (\tilde{u}_i^{n+1} u^n) = \frac{1}{|D_\sigma|} \sum_{\epsilon \in \bar{\mathcal{E}}(D_\sigma)} F_{\sigma,\epsilon} (\tilde{u}_i^{n+1})_\epsilon, $$

where $F_{\sigma,\epsilon}$ is an approximation of the mass flux through $\epsilon$ computed from the fluxes through the primal faces (i.e. for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}(K)$, the quantities $F_{K,\sigma} = \rho |\sigma| u_{K,\sigma}$ with $u_{K,\sigma}$ given by (24a) and (29) for the RT and MAC schemes respectively) in such a way that a divergence-free constraint is also satisfied on the dual cells $D_\sigma$ (i.e. $\sum_{\epsilon \in \bar{\mathcal{E}}(D_\sigma)} F_{\sigma,\epsilon} = 0$). A construction of the dual mass fluxes may be found in Appendix A of [36] (RT approximation) and Section 2 of [24] (MAC scheme). The approximation of the $i$th component of the velocity at the dual face, $(\tilde{u}_i^{n+1})_\epsilon$, may be chosen upwind or centered.

Finally, we need to ensure the duality of the velocity divergence and pressure gradient, together with the duality relations (14), (16) and (20). All these properties are obtained “by construction” when the associated terms in the momentum balance equation are derived from a variational formulation of this latter equation, as is done with the RT discretization (see Sect. 3.3.1 below). For the MAC scheme, the pressure gradient-velocity divergence duality is standard. The definition of these operators may be found, for instance, in Section 2 of [24]
**Correction step** – Solve for $p^{n+1}$, $u^{n+1}$ and $c^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \qquad \frac{\rho}{\delta t} (u_{\sigma,i}^{n+1} - \tilde{u}_{\sigma,i}^{n+1}) + \nabla_{\sigma,i} (p^{n+1} - p^n) - \text{div}_{\sigma,i} \tau_p^{n+1} + \text{div}_{\sigma,i} \tau_p^n = 0, \qquad (21b) $$

$$ \forall K \in \mathcal{M}, \quad \text{div}_K(u^{n+1}) = 0, \qquad (21c) $$

$$ \forall K \in \mathcal{M}, \qquad \frac{1}{\delta t} (c_K^{n+1} - c_K^n) + \text{div}_K(u^n c^{n+1}) = \nabla_K u^{n+1} c_K^{n+1} + (c_K^{n+1})^t (\nabla_K u^{n+1})^t - \frac{1}{\lambda} (c_K^{n+1} - I_d). \quad (21d) $$

In the prediction step and in the velocity correction step, for $k = n$ and $k = n + 1$, the tensor $\tau_p^k$ is discretized on the primal cells and is computed as a function of the conformation tensor as

$$ (\tau_p^k)_K = \frac{\eta_p}{\lambda} (c_K^k - I_d), \quad \text{for } K \in \mathcal{M}. $$

The divergence operator applied to this tensor is defined in Section 3.3.

To build a fully discrete scheme that satisfies the same stability estimate as its time-discrete counterpart, we need to consider all the properties that have been used in the proof of Proposition 3.1 and satisfy them at the discrete level by a suitable definition of the operators in the scheme. First, we require that the transport equation for the tensor $c$ preserves its symmetry and its positive-definiteness. The former property is easy to satisfy by solving a “symmetrized version” of the right-hand side (i.e. a discretization of $\nabla u c + c^t (\nabla u)^t$ instead of $\nabla u c + c (\nabla u)^t$), as is done in equation (21d). Once the equation is symmetrized, if $c^n$ is symmetric, for $1 \le i, j \le d$, the discrete equation for $c_{i,j}^{n+1}$ is strictly the same as the equation for $c_{j,i}^{n+1}$. Conversely, a symmetric tensor solution to (21d) is also a solution of the original equation and thus switching to (21d) preserves the consistency of the scheme. Moreover, it is proved in Appendix A that a first-order backward upwind scheme for (21d) preserves the positive definiteness of $c$, under a condition for the time-step depending on the $L^\infty$-norm of $\nabla u^{n+1}$. Finally, for this equation and this convection scheme, we need a discrete analogue of (19) which is proved in Appendix C.

Let us now turn to the momentum balance equation. A finite-volume discretization of the convection term for the RT discretization may be found, in the more general context of variable density flows, in Section 4.2 of [35]. The discrete kinetic energy balance (i.e. the discrete equivalent of (13)) is established in Lemma 7.2 of [35]. For the MAC scheme, a discretization of the convection term is proposed in Section 2, equation (29) of [24] and its $L^2$-stability (i.e. as before, the discrete equivalent of (13)) is shown in Lemma 3.6 of [24]. The convection operator takes the following general form

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \quad \rho \text{div}_\sigma (\tilde{u}_i^{n+1} u^n) = \frac{1}{|D_\sigma|} \sum_{\epsilon \in \bar{\mathcal{E}}(D_\sigma)} F_{\sigma, \epsilon} (\tilde{u}_i^{n+1})_\epsilon, $$

where $F_{\sigma, \epsilon}$ is an approximation of the mass flux through $\epsilon$ computed from the fluxes through the primal faces (i.e. for $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}(K)$, the quantities $F_{K, \sigma} = \rho |\sigma| u_{K, \sigma}$ with $u_{K, \sigma}$ given by (24a) and (29) for the RT and MAC schemes respectively) in such a way that a divergence-free constraint is also satisfied on the dual cells $D_\sigma$ (i.e. $\sum_{\epsilon \in \bar{\mathcal{E}}(D_\sigma)} F_{\sigma, \epsilon} = 0$). A construction of the dual mass fluxes may be found in Appendix A of [36] (RT approximation) and Section 2 of [24] (MAC scheme). The approximation of the $i$th component of the velocity at the dual face, $(\tilde{u}_i^{n+1})_\epsilon$, may be chosen upwind or centered.

Finally, we need to ensure the duality of the velocity divergence and pressure gradient, together with the duality relations (14), (16) and (20). All these properties are obtained “by construction” when the associated terms in the momentum balance equation are derived from a variational formulation of this latter equation, as is done with the RT discretization (see Sect. 3.3.1 below). For the MAC scheme, the pressure gradient-velocity divergence duality is standard. The definition of these operators may be found, for instance, in Section 2 of [24]
and a proof of their duality is given by Lemma 2.4 of [24]. For the viscous diffusion term, we propose in [25] a construction that ensures the positivity of the dissipation, i.e. a discrete analogue of Relation (14). It consists in obtaining the usual finite-volume discretization of this term by a variational technique, thanks to a discrete $H^1$-inner product defined for “finite-volume discrete functions”. For the sake of completeness, we recall these developments in Section 3.3.2 below, and give the definition of the velocity divergence, velocity diffusion, and pressure gradient operators. Then, in the same section, we extend this technique to cope with the term $\text{div}_{\sigma,i} \tau_p$, which also yields a definition of the velocity gradient used at the right-hand side of the transport equation (21d) of the conformation tensor.

### 3.3. The discrete total stress divergence term

The aim of this section is to define the divergence term $\text{div}_{\sigma,i}(\mathbf{T})$ of the total stress tensor $\mathbf{T} = -p\mathbf{I}_d + \tau_s + \tau_p$. We want this quantity to satisfy a discrete counterpart of the duality relations (14), (16) and (20) that, as shown in the previous sections, are crucial for the preservation of the free-energy estimate by the scheme. Gathering these three relations in a single one, we thus search for an identity of the form:

$$ -\sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_s^{(i)}} |D_{\sigma}| \text{div}_{\sigma,i} \mathbf{T} \, u_{\sigma,i} = a_{\mathcal{M}}(\mathbf{u}, \mathbf{u}) + \sum_{K \in \mathcal{M}} |K|(p_K \text{div}_K \mathbf{u} + \tau_p(c_K) : \nabla_K \mathbf{u}), \quad (22) $$

where $a_{\mathcal{M}}(\mathbf{u}, \mathbf{u})$ stands for non-negative mesh-dependent bilinear form associated to the viscous dissipation and the terms $\text{div}_K \mathbf{u}$ and $\nabla_K \mathbf{u}$, to be used in the discretization of equations (21c) and (21d) respectively, have to be defined. The derivation of $\text{div}_{\sigma,i} \mathbf{T}$ is different for the RT discretization and for the MAC scheme, and we deal with these two cases separately in the two following sub-sections.

#### 3.3.1. Unstructured discretization

For the RT discretization, we use the usual finite element discretization of the stress tensor divergence term

$$ -\text{div}_{\sigma,i} \mathbf{T} = \frac{1}{|D_{\sigma}|} \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \varphi_{\sigma}^{(i)} \, d\mathbf{x}, \\ \mathbf{T}(\mathbf{x}) = \eta_s \left( (\nabla \mathbf{u}(\mathbf{x}) + (\nabla \mathbf{u}(\mathbf{x}))^t)^t + p_K \mathbf{I}_d + \tau_p(c_K) \right) \quad \text{for } \mathbf{x} \in K, \quad (23) $$

where $\varphi_{\sigma}^{(i)}$ stands for the vector-valued finite element shape function associated with the $i$th component of the velocity and the face $\sigma$. By definition of the RT finite elements, this shape function is $\varphi_{\sigma}e^{(i)}$, where $\varphi_{\sigma}$ is the real-valued function of the approximation space whose mean value is 1 over $\sigma$ and 0 over the other faces of the mesh. We recall that, for the so-called parametric variant of the RT element, the discrete functional space is obtained through the usual $Q_1$-mapping from the space spanned by the set of functions $\{1, x_1, x_2, x_1^2 - x_2^2\}$ on the reference $(0,1)^2$ square ($d=2$) and $\{1, x_1, x_2, x_3, x_1^2 - x_2^2, x_2^2 - x_3^2\}$ on the reference $(0,1)^3$ cube ($d=3$).

By a simple reordering of the sums, equation (22) reads

$$ -\sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} |D_{\sigma}| \text{div}_{\sigma,i} \mathbf{T} \, u_{\sigma,i} = \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} u_{\sigma,i} \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \varphi_{\sigma}^{(i)} \, d\mathbf{x} \\ = \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \left( \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} u_{\sigma,i} \varphi_{\sigma}^{(i)} \right) \, d\mathbf{x} = \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \mathbf{u} \, d\mathbf{x} \\ = \sum_{K \in \mathcal{M}} |K| (\mathbf{T} : \nabla \mathbf{u})_K. $$
and a proof of their duality is given by Lemma 2.4 of [24]. For the viscous diffusion term, we propose in [25] a construction that ensures the positivity of the dissipation, i.e. a discrete analogue of Relation (14). It consists in obtaining the usual finite-volume discretization of this term by a variational technique, thanks to a discrete $H^1$-inner product defined for “finite-volume discrete functions”. For the sake of completeness, we recall these developments in Section 3.3.2 below, and give the definition of the velocity divergence, velocity diffusion, and pressure gradient operators. Then, in the same section, we extend this technique to cope with the term $\text{div}_{\sigma,i} \tau_p$, which also yields a definition of the velocity gradient used at the right-hand side of the transport equation (21d) of the conformation tensor.

### 3.3. The discrete total stress divergence term

The aim of this section is to define the divergence term $\text{div}_{\sigma,i}(\mathbf{T})$ of the total stress tensor $\mathbf{T} = -p\mathbf{I}_d + \tau_s + \tau_p$. We want this quantity to satisfy a discrete counterpart of the duality relations (14), (16) and (20) that, as shown in the previous sections, are crucial for the preservation of the free-energy estimate by the scheme. Gathering these three relations in a single one, we thus search for an identity of the form:

$$ -\sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_s^{(i)}} |D_{\sigma}| \text{div}_{\sigma,i} \mathbf{T} \, u_{\sigma,i} = a_{\mathcal{M}}(\mathbf{u}, \mathbf{u}) + \sum_{K \in \mathcal{M}} |K|(p_K \text{div}_K \mathbf{u} + \tau_p(c_K) : \nabla_K \mathbf{u}), \quad (22) $$

where $a_{\mathcal{M}}(\mathbf{u}, \mathbf{u})$ stands for non-negative mesh-dependent bilinear form associated to the viscous dissipation and the terms $\text{div}_K \mathbf{u}$ and $\nabla_K \mathbf{u}$, to be used in the discretization of equations (21c) and (21d) respectively, have to be defined. The derivation of $\text{div}_{\sigma,i} \mathbf{T}$ is different for the RT discretization and for the MAC scheme, and we deal with these two cases separately in the two following sub-sections.

#### 3.3.1. Unstructured discretization

For the RT discretization, we use the usual finite element discretization of the stress tensor divergence term

$$ -\text{div}_{\sigma,i} \mathbf{T} = \frac{1}{|D_{\sigma}|} \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \varphi_{\sigma}^{(i)} \, d\mathbf{x}, \\ \mathbf{T}(\mathbf{x}) = \eta_s \left( (\nabla \mathbf{u}(\mathbf{x}) + (\nabla \mathbf{u}(\mathbf{x}))^t)^t + p_K \mathbf{I}_d + \tau_p(c_K) \right) \quad \text{for } \mathbf{x} \in K, \quad (23) $$

where $\varphi_{\sigma}^{(i)}$ stands for the vector-valued finite element shape function associated with the $i$th component of the velocity and the face $\sigma$. By definition of the RT finite elements, this shape function is $\varphi_{\sigma} e^{(i)}$, where $\varphi_{\sigma}$ is the real-valued function of the approximation space whose mean value is 1 over $\sigma$ and 0 over the other faces of the mesh. We recall that, for the so-called parametric variant of the RT element, the discrete functional space is obtained through the usual $Q_1$-mapping from the space spanned by the set of functions $\{1, x_1, x_2, x_1^2 - x_2^2\}$ on the reference $(0,1)^2$ square ($d=2$) and $\{1, x_1, x_2, x_3, x_1^2 - x_2^2, x_2^2 - x_3^2\}$ on the reference $(0,1)^3$ cube ($d=3$).

By a simple reordering of the sums, equation (22) reads

$$ -\sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} |D_{\sigma}| \text{div}_{\sigma,i} \mathbf{T} \, u_{\sigma,i} = \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} u_{\sigma,i} \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \varphi_{\sigma}^{(i)} \, d\mathbf{x} \\ = \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \left( \sum_{i=1}^{d} \sum_{\sigma \in \mathcal{E}_{\text{int}}} u_{\sigma,i} \varphi_{\sigma}^{(i)} \right) \, d\mathbf{x} = \sum_{K \in \mathcal{M}} \int_K \mathbf{T} : \nabla \mathbf{u} \, d\mathbf{x} \\ = \sum_{K \in \mathcal{M}} |K| (\mathbf{T} : \nabla \mathbf{u})_K. $$
Considering separately the pressure gradient, the diffusion and the elastic extra-stress terms in $\mathcal{T}$, we obtain the following expression for the discrete operators featured in the scheme (21), dropping the time index for short

$$ \forall K \in \mathcal{M}, \qquad \operatorname{div}_K \mathbf{u} = \frac{1}{|K|} \int_K \operatorname{div}(\mathbf{u}) \, d\mathbf{x} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K)} |\sigma| u_{K,\sigma}, \qquad (24a) $$
with $u_{K,\sigma} = \mathbf{u}_\sigma \cdot \mathbf{n}_{K,\sigma}$,

$$ \forall \sigma \in \mathcal{E}_{\text{int}}, \sigma = K|L, \qquad \nabla_{\sigma} p = \frac{1}{|D_{\sigma}|} \sum_{i=1}^{d} \sum_{M \in \mathcal{M}} p_{M} \int_{M} \nabla \varphi_{\sigma}^{(i)} d\mathbf{x} = \frac{|\sigma|}{|D_{\sigma}|} (p_{L} - p_{K}) \mathbf{n}_{K,\sigma}, \qquad (24b) $$

$$ \text{For } 1 \le i \le d, \forall \sigma \in \mathcal{E}_{\text{int}}, \qquad \operatorname{div}_{\sigma,i} \tau_s(\tilde{u}) = \frac{1}{|D_{\sigma}|} \int_{\Omega} \tau_s(\tilde{u}) : \nabla \varphi_{\sigma}^{(i)} d\mathbf{x}, \qquad (24c) $$

$$ \text{For } 1 \le i \le d, \forall \sigma \in \mathcal{E}_{\text{int}}, \qquad \operatorname{div}_{\sigma,i} \tau_p = \frac{1}{|D_{\sigma}|} \sum_{K \in \mathcal{M}} (\tau_p)_K : \int_K \nabla \varphi_{\sigma}^{(i)} d\mathbf{x}, \qquad (24d) $$

$$ \forall K \in \mathcal{M}, \qquad \nabla_K \mathbf{u} = \frac{1}{|K|} \int_K \nabla \mathbf{u} \, d\mathbf{x}. \qquad (24e) $$

### 3.3.2. MAC scheme

For the MAC scheme, the strategy implemented here is to mimick the RT case, i.e. to recast the momentum balance equation under a weak form. We define a discrete gradient operator acting on the velocity discrete functions (i.e. piecewise-constant functions over the dual cells). Each component of this discrete gradient is itself a piecewise-constant function, based on a new partition of the computational domain. In two-dimensions, the discrete analogues of $\partial_x u_x$ and $\partial_y u_y$ are defined over the primal cells, while the other components of the velocity gradients are defined over volumes associated with the mesh vertices. With this gradient operator defined, the machinery of standard Galerkin techniques applies for the discretization of the momentum balance. The discrete velocity function space is shown to be spanned by a set of shape functions, and each discrete equation is derived by using these shape functions as test functions. Thanks to this construction, the scheme discrete operators natively satisfy the duality Relation (22) (they are indeed obtained from this relation itself), and have formally the same expression than in the finite element case, i.e. obey formula similar to System (24) above.

This approach was followed in [25] for Newtonian fluids. We briefly recall here these results and use the same technique to cope with the divergence of the elastic extra-stress tensor. For short, we only address here the two-dimensional case and we refer to [25] for the (natural) extension to the three-dimensional case.

**The discrete velocity gradient.** The developments of this section rely on the geometric properties of structured meshes and we need to momentarily switch to the notations that are specific to this context. We thus suppose that the primal mesh reads $\mathcal{M} = \{K_{i,j}, (i,j) \in \mathcal{J}\}$, with $\mathcal{J} \subset \mathbb{N}^2$ with notations for geometrical quantities that are given in Figures 3–5. We detail here the discretization of the terms associated with the first component of the velocity, the adaptation to the second component being straightforward. In the interior of the computational domain, the discrete partial derivatives of this velocity component are defined as follows.

- For $(i,j) \in \mathcal{J}$, let the primal cells be $K_{i,j} = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}) \times (y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}})$. The discrete derivative involved in the divergence (so, for the velocity $x$-component, only $\partial_x^M u^x$) is defined over the primal cell by, $\forall (i,j) \in \mathcal{J}$ and $\forall \mathbf{x} \in K_{i,j}$,

$$ \partial_x^M u^x(\mathbf{x}) = \frac{u_{i+\frac{1}{2},j}^x - u_{i-\frac{1}{2},j}^x}{h_i^x}. \qquad (25) $$

- For the other derivatives (so, for the velocity $x$-component, only $\partial_y^M u^x$), we introduce a fourth mesh which is vertex-centred, and we denote by $K^{xy}$ the generic cell of this new mesh, with $K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy} = (x_{i-1}, x_i) \times$
Considering separately the pressure gradient, the diffusion and the elastic extra-stress terms in $\mathcal{T}$, we obtain the following expression for the discrete operators featured in the scheme (21), dropping the time index for short

$$ \forall K \in \mathcal{M}, \qquad \operatorname{div}_K \mathbf{u} = \frac{1}{|K|} \int_K \operatorname{div}(\mathbf{u}) \, d\mathbf{x} = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K)} |\sigma| u_{K,\sigma}, \qquad (24a) $$
with $u_{K,\sigma} = \mathbf{u}_\sigma \cdot \mathbf{n}_{K,\sigma}$,

$$ \forall \sigma \in \mathcal{E}_{\text{int}}, \sigma = K|L, \qquad \nabla_{\sigma} p = \frac{1}{|D_{\sigma}|} \sum_{i=1}^{d} \sum_{M \in \mathcal{M}} p_{M} \int_{M} \nabla \varphi_{\sigma}^{(i)} d\mathbf{x} = \frac{|\sigma|}{|D_{\sigma}|} (p_{L} - p_{K}) \mathbf{n}_{K,\sigma}, \qquad (24b) $$

$$ \text{For } 1 \le i \le d, \forall \sigma \in \mathcal{E}_{\text{int}}, \qquad \operatorname{div}_{\sigma,i} \tau_s(\tilde{u}) = \frac{1}{|D_{\sigma}|} \int_{\Omega} \tau_s(\tilde{u}) : \nabla \varphi_{\sigma}^{(i)} d\mathbf{x}, \qquad (24c) $$

$$ \text{For } 1 \le i \le d, \forall \sigma \in \mathcal{E}_{\text{int}}, \qquad \operatorname{div}_{\sigma,i} \tau_p = \frac{1}{|D_{\sigma}|} \sum_{K \in \mathcal{M}} (\tau_p)_K : \int_K \nabla \varphi_{\sigma}^{(i)} d\mathbf{x}, \qquad (24d) $$

$$ \forall K \in \mathcal{M}, \qquad \nabla_K \mathbf{u} = \frac{1}{|K|} \int_K \nabla \mathbf{u} \, d\mathbf{x}. \qquad (24e) $$

### 3.3.2. MAC scheme

For the MAC scheme, the strategy implemented here is to mimick the RT case, i.e. to recast the momentum balance equation under a weak form. We define a discrete gradient operator acting on the velocity discrete functions (i.e. piecewise-constant functions over the dual cells). Each component of this discrete gradient is itself a piecewise-constant function, based on a new partition of the computational domain. In two-dimensions, the discrete analogues of $\partial_x u_x$ and $\partial_y u_y$ are defined over the primal cells, while the other components of the velocity gradients are defined over volumes associated with the mesh vertices. With this gradient operator defined, the machinery of standard Galerkin techniques applies for the discretization of the momentum balance. The discrete velocity function space is shown to be spanned by a set of shape functions, and each discrete equation is derived by using these shape functions as test functions. Thanks to this construction, the scheme discrete operators natively satisfy the duality Relation (22) (they are indeed obtained from this relation itself), and have formally the same expression than in the finite element case, i.e. obey formula similar to System (24) above.

This approach was followed in [25] for Newtonian fluids. We briefly recall here these results and use the same technique to cope with the divergence of the elastic extra-stress tensor. For short, we only address here the two-dimensional case and we refer to [25] for the (natural) extension to the three-dimensional case.

**The discrete velocity gradient.** The developments of this section rely on the geometric properties of structured meshes and we need to momentarily switch to the notations that are specific to this context. We thus suppose that the primal mesh reads $\mathcal{M} = \{K_{i,j}, (i,j) \in \mathcal{J}\}$, with $\mathcal{J} \subset \mathbb{N}^2$ with notations for geometrical quantities that are given in Figures 3–5. We detail here the discretization of the terms associated with the first component of the velocity, the adaptation to the second component being straightforward. In the interior of the computational domain, the discrete partial derivatives of this velocity component are defined as follows.

- For $(i,j) \in \mathcal{J}$, let the primal cells be $K_{i,j} = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}) \times (y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}})$. The discrete derivative involved in the divergence (so, for the velocity $x$-component, only $\partial_x^M u^x$) is defined over the primal cell by, $\forall (i,j) \in \mathcal{J}$ and $\forall \mathbf{x} \in K_{i,j}$,

$$ \partial_x^M u^x(\mathbf{x}) = \frac{u_{i+\frac{1}{2},j}^x - u_{i-\frac{1}{2},j}^x}{h_i^x}. \qquad (25) $$

- For the other derivatives (so, for the velocity $x$-component, only $\partial_y^M u^x$), we introduce a fourth mesh which is vertex-centred, and we denote by $K^{xy}$ the generic cell of this new mesh, with $K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy} = (x_{i-1}, x_i) \times$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Unknown and dual cell for the *x*-component of the velocity, notations for staggered discretizations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Discrete partial derivatives of the *x*-component of the velocity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Boundary conditions, *x*-component of the velocity.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Unknown and dual cell for the *x*-component of the velocity, notations for staggered discretizations.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4. Discrete partial derivatives of the *x*-component of the velocity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Boundary conditions, *x*-component of the velocity.</figcaption></figure>
$(y_{j-1}, y_j)$. We define $\mathcal{J}^{xy}$ the set of pairs $(i, j)$ such that $x_{i-\frac{1}{2},j-\frac{1}{2}}$ is an internal vertex of the mesh (i.e. a vertex which does not lie on the boundary). Then, for $(i, j) \in \mathcal{J}^{xy}$ and $\forall x \in K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy}$:

$$ \partial_y^M u^x(\mathbf{x}) = \frac{u_{i-\frac{1}{2},j}^x - u_{i-\frac{1}{2},j-1}^x}{h_{j-\frac{1}{2}}^y}. \qquad (26) $$

The only necessary extension of this definition to cope with boundaries concerns the definition of $\partial_y^M u^x$ over a half vertex-centered cell associated with a vertex lying on a horizontal boundary (see Fig. 5). In this case, we use the usual “fictitious cell trick” in order to apply Relation (26): an external cell, of zero y-dimension, is added to the mesh and the horizontal velocity in this cell is set to the prescribed Dirichlet value, or to zero for the test functions defined below. Extending these definitions to the y-component of the velocity, the discrete diffusion tensor can be defined as

$$ \nabla^M \mathbf{u} = \begin{bmatrix} \partial_x^M u^x & \partial_y^M u^x \\ \partial_x^M u^y & \partial_y^M u^y \end{bmatrix}, \quad \tau_s^M(\mathbf{u}) = \eta_s \left( \nabla^M \mathbf{u} + (\nabla^M \mathbf{u})^t \right). \qquad (27) $$

**Finite-volume test functions.** Let us denote by $\mathcal{J}^x \subset \mathbb{N}^2$ (resp. $\mathcal{J}^y \subset \mathbb{N}^2$) the set of pairs $(i, j)$ such that $x_{i-\frac{1}{2},j}$ (resp. $x_{i,j-\frac{1}{2}}$) is the mass center of a vertical (resp. horizontal) internal face of the mesh. For $(i, j) \in \mathcal{J}^x$, we denote by $\phi^{x,(i-\frac{1}{2},j)}$ the test function associated with the degree of freedom of the x-component of the velocity located at $x_{i-\frac{1}{2},j}$. This discrete function is defined as

$$ (\phi^{x,(i-\frac{1}{2},j)})_{k-\frac{1}{2},\ell}^x = \delta_k^i \delta_\ell^j, \quad \forall (k, \ell) \in \mathcal{J}^x \text{ and } (\phi^{x,(i-\frac{1}{2},j)})_{k,\ell-\frac{1}{2}}^y = 0, \quad \forall (k, \ell) \in \mathcal{J}^y. $$

Its non-zero partial derivatives are $\partial_x^M \phi^{x,(i-\frac{1}{2},j)}$ and $\partial_y^M \phi^{x,(i-\frac{1}{2},j)}$ with

$$ \partial_x^M \phi^{x,(i-\frac{1}{2},j)} = \begin{cases} \frac{1}{h_{i-1}^x} & \text{over } K_{i-1,j}, \\ -\frac{1}{h_i^x} & \text{over } K_{i,j}, \\ 0 & \text{elsewhere,} \end{cases} \qquad \partial_y^M \phi^{x,(i-\frac{1}{2},j)} = \begin{cases} \frac{1}{h_{j-\frac{1}{2}}^y} & \text{over } K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy}, \\ -\frac{1}{h_{j+\frac{1}{2}}^y} & \text{over } K_{i-\frac{1}{2},j+\frac{1}{2}}^{xy}, \\ 0 & \text{elsewhere.} \end{cases} $$

Since the velocity is prescribed on the boundary, no equation is written on the half-dual cells associated with external faces, so no definition is required for the corresponding test functions.

**Discrete pressure gradient, velocity divergence and viscous diffusion.** Let us begin with the pressure gradient. Identifying $p$ with its associated piecewise constant function, we have

$$ \forall (i, j) \in \mathcal{J}^x, \quad (\nabla p)_{i-\frac{1}{2},j}^x = \frac{-1}{|K_{i-\frac{1}{2},j}|} \int_{\Omega} p \partial_x^M \phi^{x,(i-\frac{1}{2},j)} \, d\mathbf{x}. \qquad (28) $$

Using the expression of the discrete derivatives of the test functions, we get

$$ (\nabla p)_{i-\frac{1}{2},j}^x = \frac{-1}{|K_{i-\frac{1}{2},j}|} \left( (h_{i-1}^x h_j^y) p_{i-1,j} + \frac{1}{h_{i-1}^x} + (h_i^x h_j^y) p_{i,j} - \frac{1}{h_i^x} \right) = \frac{p_{i,j} - p_{i-1,j}}{h_{i-\frac{1}{2}}^x}, $$

and we recover the usual differential quotient for the pressure gradient of the MAC scheme. Note that, returning to the general notations used in this paper, this expression may also be recast as

$$ \text{for } 1 \le i \le d, \forall \sigma \in \mathcal{E}_S^{(i)}, \sigma = K|L, \quad \nabla_{\sigma,i} p = \frac{|\sigma|}{|D_\sigma|} (p_L - p_K) (n_{K,\sigma})_i, $$
$(y_{j-1}, y_j)$. We define $\mathcal{J}^{xy}$ the set of pairs $(i, j)$ such that $x_{i-\frac{1}{2},j-\frac{1}{2}}$ is an internal vertex of the mesh (i.e. a vertex which does not lie on the boundary). Then, for $(i, j) \in \mathcal{J}^{xy}$ and $\forall x \in K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy}$:

$$ \partial_y^M u^x(\mathbf{x}) = \frac{u_{i-\frac{1}{2},j}^x - u_{i-\frac{1}{2},j-1}^x}{h_{j-\frac{1}{2}}^y}. \qquad (26) $$

The only necessary extension of this definition to cope with boundaries concerns the definition of $\partial_y^M u^x$ over a half vertex-centered cell associated with a vertex lying on a horizontal boundary (see Fig. 5). In this case, we use the usual “fictitious cell trick” in order to apply Relation (26): an external cell, of zero y-dimension, is added to the mesh and the horizontal velocity in this cell is set to the prescribed Dirichlet value, or to zero for the test functions defined below. Extending these definitions to the y-component of the velocity, the discrete diffusion tensor can be defined as

$$ \nabla^M \mathbf{u} = \begin{bmatrix} \partial_x^M u^x & \partial_y^M u^x \\ \partial_x^M u^y & \partial_y^M u^y \end{bmatrix}, \quad \tau_s^M(\mathbf{u}) = \eta_s \left( \nabla^M \mathbf{u} + (\nabla^M \mathbf{u})^t \right). \qquad (27) $$

**Finite-volume test functions.** Let us denote by $\mathcal{J}^x \subset \mathbb{N}^2$ (resp. $\mathcal{J}^y \subset \mathbb{N}^2$) the set of pairs $(i, j)$ such that $x_{i-\frac{1}{2},j}$ (resp. $x_{i,j-\frac{1}{2}}$) is the mass center of a vertical (resp. horizontal) internal face of the mesh. For $(i, j) \in \mathcal{J}^x$, we denote by $\phi^{x,(i-\frac{1}{2},j)}$ the test function associated with the degree of freedom of the x-component of the velocity located at $x_{i-\frac{1}{2},j}$. This discrete function is defined as

$$ (\phi^{x,(i-\frac{1}{2},j)})_{k-\frac{1}{2},\ell}^x = \delta_k^i \delta_\ell^j, \quad \forall (k, \ell) \in \mathcal{J}^x \text{ and } (\phi^{x,(i-\frac{1}{2},j)})_{k,\ell-\frac{1}{2}}^y = 0, \quad \forall (k, \ell) \in \mathcal{J}^y. $$

Its non-zero partial derivatives are $\partial_x^M \phi^{x,(i-\frac{1}{2},j)}$ and $\partial_y^M \phi^{x,(i-\frac{1}{2},j)}$ with

$$ \partial_x^M \phi^{x,(i-\frac{1}{2},j)} = \begin{cases} \frac{1}{h_{i-1}^x} & \text{over } K_{i-1,j}, \\ -\frac{1}{h_i^x} & \text{over } K_{i,j}, \\ 0 & \text{elsewhere,} \end{cases} \qquad \partial_y^M \phi^{x,(i-\frac{1}{2},j)} = \begin{cases} \frac{1}{h_{j-\frac{1}{2}}^y} & \text{over } K_{i-\frac{1}{2},j-\frac{1}{2}}^{xy}, \\ -\frac{1}{h_{j+\frac{1}{2}}^y} & \text{over } K_{i-\frac{1}{2},j+\frac{1}{2}}^{xy}, \\ 0 & \text{elsewhere.} \end{cases} $$

Since the velocity is prescribed on the boundary, no equation is written on the half-dual cells associated with external faces, so no definition is required for the corresponding test functions.

**Discrete pressure gradient, velocity divergence and viscous diffusion.** Let us begin with the pressure gradient. Identifying $p$ with its associated piecewise constant function, we have

$$ \forall (i, j) \in \mathcal{J}^x, \quad (\nabla p)_{i-\frac{1}{2},j}^x = \frac{-1}{|K_{i-\frac{1}{2},j}|} \int_{\Omega} p \partial_x^M \phi^{x,(i-\frac{1}{2},j)} \mathrm{d}\mathbf{x}. \qquad (28) $$

Using the expression of the discrete derivatives of the test functions, we get

$$ (\nabla p)_{i-\frac{1}{2},j}^x = \frac{-1}{|K_{i-\frac{1}{2},j}|} \left( (h_{i-1}^x h_j^y) p_{i-1,j} + \frac{1}{h_{i-1}^x} (h_i^x h_j^y) p_{i,j} \frac{-1}{h_i^x} \right) = \frac{p_{i,j} - p_{i-1,j}}{h_{i-\frac{1}{2}}^x}, $$

and we recover the usual differential quotient for the pressure gradient of the MAC scheme. Note that, returning to the general notations used in this paper, this expression may also be recast as

$$ \text{for } 1 \le i \le d, \forall \sigma \in \mathcal{E}_S^{(i)}, \sigma = K|L, \nabla_{\sigma,i} p = \frac{|\sigma|}{|D_\sigma|} (p_L - p_K) (n_{K,\sigma})_i, $$
which is a definition similar to the RT discretization, equation (24b). The velocity divergence operator immediately follows from the definition of the velocity discrete partial derivatives

$$ \forall K \in \mathcal{M}, \quad \mathrm{div}_K(\mathbf{u}) = \frac{1}{|K|} \int_K \mathrm{div}(\mathbf{u}) \, d\mathbf{x}, $$

so, since the partial derivatives of the velocity involved in the divergence are constant over each primal cell,

$$ \forall (i,j) \in \mathcal{J}, \quad \mathrm{div}_{K_{i,j}}(\mathbf{u}) = \frac{u_{i+\frac{1}{2},j}^x - u_{i-\frac{1}{2},j}^x}{h_i^x} + \frac{u_{i,j+\frac{1}{2}}^y - u_{i,j-\frac{1}{2}}^y}{h_j^y}, $$

which is, once again, the standard definition for the MAC scheme. Note that, to get similar expressions for the RT and MAC discretizations, we may also recast this equation as

$$ \forall K \in \mathcal{M}, \quad \mathrm{div}_K(\mathbf{u}) = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K)} |\sigma| u_{K,\sigma}, \quad \text{with, for } \sigma \in \mathcal{E}^{(i)}, u_{K,\sigma} = u_{\sigma} n_{K,\sigma} \cdot e^{(i)}. \quad (29) $$

Finally, the discrete divergence of the viscous stress tensor is defined by the following weak formulation

$$ \forall (i,j) \in \mathcal{J}^x, \quad -(\mathrm{div}\tau_s(\tilde{\mathbf{u}}))_{i-\frac{1}{2},j}^x = \frac{1}{|K_{i-\frac{1}{2},j}^x|} \int_{\Omega} \tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) : \nabla^{\mathcal{M}} \phi^{x,(i-\frac{1}{2},j)} d\mathbf{x}, \quad (30) $$

By a standard reordering of the sums, we get

$$ \sum_{i=1}^{d} \sum_{\mathbf{x} \in \mathcal{E}_g^{(i)}} |D_{\sigma}| \tilde{u}_{\sigma,i} \mathrm{div}_{\sigma,i} \tau_s(\tilde{\mathbf{u}}) = \int_{\Omega} \tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) : \nabla^{\mathcal{M}} \tilde{\mathbf{u}} \, d\mathbf{x}, \quad (31) $$

An important consequence of this relation, with the definition of $\tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) = \eta_s(\nabla^{\mathcal{M}}\tilde{\mathbf{u}} + (\nabla^{\mathcal{M}}\tilde{\mathbf{u}})^t)$, is that, thanks to an elementary property of the tensors contraction, we may replace the second tensor $\nabla^{\mathcal{M}}\tilde{\mathbf{u}}$ in (31) by its symmetrization $(\nabla^{\mathcal{M}}\tilde{\mathbf{u}} + (\nabla^{\mathcal{M}}\tilde{\mathbf{u}})^t)/2$. We thus obtain that the dissipation (i.e. the integral in (31)) is non-negative, which is a discrete counterpart of (14). In addition, the duality Relation (20) requires the discretization of $\nabla_K \mathbf{u}$ in the right-hand side of the transport equation for the conformation tensor (21d) to be

$$ \nabla_K \mathbf{u} = \frac{1}{|K|} \int_K \nabla^{\mathcal{M}} \mathbf{u} \, d\mathbf{x}. $$

This weak formulation of the diffusion term (31) is shown in [25] to lead to the usual diffusion operator when the viscosity does not depend on space. The same is true (i.e. the weak formulation yields the usual conservative finite-volume scheme) for a space-dependent viscosity only when this quantity is considered as a piecewise-constant function over the same control volumes as the discrete velocity gradient components. For instance, for two-dimensional spaces, in the compression terms $(\tau_s^{\mathcal{M}})_{x,x}$ and $(\tau_s^{\mathcal{M}})_{y,y}$, the viscosity must be associated with the primal cells, while, for the shear components $(\tau_s^{\mathcal{M}})_{x,y}$ and $(\tau_s^{\mathcal{M}})_{y,x}$, it must be associated with the vertices (which, when constructing the scheme by the usual finite volume approach, is in fact rather natural).

**Polymeric stress tensor divergence.** Let us now derive the discretization of the divergence of the polymeric stress tensor. To this end, we first associate the discrete polymeric stress $\tau_p$ to a piecewise function over the primary cells by

$$ \forall (i,j) \in \mathcal{J}, \forall \mathbf{x} \in K_{i,j}, \quad \tau_p(\mathbf{x}) = \frac{\eta_p}{\lambda} f(c_{i,j})(c_{i,j} - \mathbf{I}_d). $$

Then we set

$$ \forall (i,j) \in \mathcal{J}^x, \quad -(\mathrm{div}\tau_p)_{i-\frac{1}{2},j}^x = \frac{1}{|K_{i-\frac{1}{2},j}^x|} \int_{\Omega} \tau_p : \nabla^{\mathcal{M}} \phi^{x,(i-\frac{1}{2},j)} \, d\mathbf{x}. $$
which is a definition similar to the RT discretization, equation (24b). The velocity divergence operator immediately follows from the definition of the velocity discrete partial derivatives

$$ \forall K \in \mathcal{M}, \quad \mathrm{div}_K(\mathbf{u}) = \frac{1}{|K|} \int_K \mathrm{div}(\mathbf{u}) \, d\mathbf{x}, $$

so, since the partial derivatives of the velocity involved in the divergence are constant over each primal cell,

$$ \forall (i,j) \in \mathcal{J}, \quad \mathrm{div}_{K_{i,j}}(\mathbf{u}) = \frac{u_{i+\frac{1}{2},j}^x - u_{i-\frac{1}{2},j}^x}{h_i^x} + \frac{u_{i,j+\frac{1}{2}}^y - u_{i,j-\frac{1}{2}}^y}{h_j^y}, $$

which is, once again, the standard definition for the MAC scheme. Note that, to get similar expressions for the RT and MAC discretizations, we may also recast this equation as

$$ \forall K \in \mathcal{M}, \quad \mathrm{div}_K(\mathbf{u}) = \frac{1}{|K|} \sum_{\sigma \in \mathcal{E}(K)} |\sigma| u_{K,\sigma}, \quad \text{with, for } \sigma \in \mathcal{E}^{(i)}, u_{K,\sigma} = u_{\sigma} n_{K,\sigma} \cdot e^{(i)}. \quad (29) $$

Finally, the discrete divergence of the viscous stress tensor is defined by the following weak formulation

$$ \forall (i,j) \in \mathcal{J}^x, \quad -(\mathrm{div}\tau_s(\tilde{\mathbf{u}}))_{i-\frac{1}{2},j}^x = \frac{1}{|K_{i-\frac{1}{2},j}^x|} \int_{\Omega} \tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) : \nabla^{\mathcal{M}} \phi^{x,(i-\frac{1}{2},j)} d\mathbf{x}, \quad (30) $$

By a standard reordering of the sums, we get

$$ \sum_{i=1}^{d} \sum_{\mathbf{x} \in \mathcal{E}_g^{(i)}} |D_{\sigma}| \tilde{u}_{\sigma,i} \mathrm{div}_{\sigma,i} \tau_s(\tilde{\mathbf{u}}) = \int_{\Omega} \tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) : \nabla^{\mathcal{M}} \tilde{\mathbf{u}} \, d\mathbf{x}, \quad (31) $$

An important consequence of this relation, with the definition of $\tau_s^{\mathcal{M}}(\tilde{\mathbf{u}}) = \eta_s(\nabla^{\mathcal{M}}\tilde{\mathbf{u}} + (\nabla^{\mathcal{M}}\tilde{\mathbf{u}})^t)$, is that, thanks to an elementary property of the tensors contraction, we may replace the second tensor $\nabla^{\mathcal{M}}\tilde{\mathbf{u}}$ in (31) by its symmetrization $(\nabla^{\mathcal{M}}\tilde{\mathbf{u}} + (\nabla^{\mathcal{M}}\tilde{\mathbf{u}})^t)/2$. We thus obtain that the dissipation (i.e. the integral in (31)) is non-negative, which is a discrete counterpart of (14). In addition, the duality Relation (20) requires the discretization of $\nabla_K \mathbf{u}$ in the right-hand side of the transport equation for the conformation tensor (21d) to be

$$ \nabla_K \mathbf{u} = \frac{1}{|K|} \int_K \nabla^{\mathcal{M}} \mathbf{u} \, d\mathbf{x}. $$

This weak formulation of the diffusion term (31) is shown in [25] to lead to the usual diffusion operator when the viscosity does not depend on space. The same is true (i.e. the weak formulation yields the usual conservative finite-volume scheme) for a space-dependent viscosity only when this quantity is considered as a piecewise-constant function over the same control volumes as the discrete velocity gradient components. For instance, for two-dimensional spaces, in the compression terms $(\tau_s^{\mathcal{M}})_{x,x}$ and $(\tau_s^{\mathcal{M}})_{y,y}$, the viscosity must be associated with the primal cells, while, for the shear components $(\tau_s^{\mathcal{M}})_{x,y}$ and $(\tau_s^{\mathcal{M}})_{y,x}$, it must be associated with the vertices (which, when constructing the scheme by the usual finite volume approach, is in fact rather natural).

**Polymeric stress tensor divergence.** Let us now derive the discretization of the divergence of the polymeric stress tensor. To this end, we first associate the discrete polymeric stress $\tau_p$ to a piecewise function over the primary cells by

$$ \forall (i,j) \in \mathcal{J}, \quad \forall \mathbf{x} \in K_{i,j}, \quad \tau_p(\mathbf{x}) = \frac{\eta_p}{\lambda} f(c_{i,j})(c_{i,j} - \mathbf{I}_d). $$

Then we set

$$ \forall (i,j) \in \mathcal{J}^x, \quad -(\mathrm{div}\tau_p)_{i-\frac{1}{2},j}^x = \frac{1}{|K_{i-\frac{1}{2},j}^x|} \int_{\Omega} \tau_p : \nabla^{\mathcal{M}} \phi^{x,(i-\frac{1}{2},j)}. $$
An easy calculation shows that this relation may be recast as a finite volume formulation, in the sense that the right-hand side may be seen as a sum over the faces of $K_{i-\frac{1}{2},j}^x$ of a discretization of the flux associated with $\text{div}\tau_p$, i.e. the integral of the first component of $\tau_p n_{K,\sigma}$ (see Fig. 3).

- On the vertical face $\varepsilon = \varepsilon_{i,j}^x$, the flux is given by

$$\int_{\varepsilon} (\tau_p n_{\varepsilon})^x = \frac{1}{h_j^y} (\tau_p)_{i,j}^{xx} = |\varepsilon_{i,j}^x| (\tau_p)_{i,j}^{xx}.$$

- The flux through the horizontal half face $\varepsilon = \varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x \cap \bar{K}_{i,j}$ reads

$$\int_{\varepsilon} (\tau_p n_{\varepsilon})^x = \frac{1}{h_{j+\frac{1}{2}}^y} \left( \frac{|K_{i,j}|}{4} (\tau_p)_{i,j}^{xy} + \frac{|K_{i,j+1}|}{4} (\tau_p)_{i,j+1}^{xy} \right) = |\varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x| \left( \alpha (\tau_p)_{i,j}^{xy} + (1-\alpha) (\tau_p)_{i,j+1}^{xy} \right),$$

with

$$\alpha = \frac{h_j^y}{2h_{j+\frac{1}{2}}^y}, \quad \text{and} \quad 1 - \alpha = \frac{h_{j+1}^y}{2h_{j+\frac{1}{2}}^y}.$$

Hence, as usual when such a duality technique is used, the approximation of the tensor at the horizontal faces may seem strange. It is a convex combination of the unknowns in the two neighbouring cells, but with coefficients that are not those which would be given by a linear interpolation. Indeed, on $\varepsilon = \varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x$, we have

$$(\tau_p)_{\varepsilon}^{xy} = \alpha (\tau_p)_{i,j}^{xy} + (1-\alpha) (\tau_p)_{i,j+1}^{xy},$$

whereas, with a linear interpolation, we would have

$$(\tau_p)_{\varepsilon}^{xy} = (1-\alpha) (\tau_p)_{i,j}^{xy} + \alpha (\tau_p)_{i,j+1}^{xy},$$

so that the formulae differ on a non-uniform mesh.

Up to now, we have only invoked stability arguments to justify the derivation of the discrete divergence of the elastic extra-stress tensor through a weak formulation. However, thanks to the consistency properties of the discrete gradient $\nabla^M$, this technique also ensures by construction the weak consistency of this discretization in the Lax-Wendroff sense, as stated in the following result for the Oldroyd-B model. For short, it is written for functions depending only on space but trivially extends to time-dependent functions.

**Lemma 3.2.** Let $(\mathcal{M}^{(m)})_{m \in \mathbb{N}}$ be sequence of structured discretizations of $\Omega$, with a space step tending to zero. Let $(c^{(m)})_{m \in \mathbb{N}}$ be a sequence of discrete conformation tensors, with $c^{(m)}$ associated with $\mathcal{M}^{(m)}$. We suppose that $(c^{(m)})_{m \in \mathbb{N}}$ weakly converges in $L^1(\Omega)^{d \times d}$ to $c$. Let $\varphi \in C_c^\infty(\Omega)^d$, and let us denote by $\varphi^{(m)}$ the interpolate of $\varphi$ in the space of velocity discrete functions defined by $\varphi^{(m)}(x) = \varphi_i(x_\sigma)$ for $x \in D_\sigma$, $\sigma \in (\varepsilon^{(m)})^{(i)}$, $1 \le i \le d$, with $x_\sigma$ the mass center of $\sigma$. Finally, for a discrete tensor $T$, i.e. a tensor-valued piecewise constant function over the primal cells of $\mathcal{M}^{(m)}$, we denote by $\text{div}_{\mathcal{M}^{(m)}} T: \Omega \to \mathbb{R}^d$ the vector-valued function defined by $(\text{div}_{\mathcal{M}^{(m)}} T)_i = \text{div}_{\sigma,i} T$ for $x \in D_\sigma$, $\sigma \in (\varepsilon^{(m)})^{(i)}$, $1 \le i \le d$. Then

$$\lim_{m \to +\infty} \int_{\Omega} \text{div}_{\mathcal{M}^{(m)}} \left( \frac{\eta_p}{\lambda} (c^{(m)} - I_d) \right) \cdot \varphi_i^{(m)} \,dx = - \int_{\Omega} \frac{\eta_p}{\lambda} (c - I_d) : \nabla \varphi \,dx.$$

*Proof.* By construction, we have

$$\int_{\Omega} \text{div}_{\mathcal{M}^{(m)}} \left( \frac{\eta_p}{\lambda} (c^{(m)} - I_d) \right) \cdot \varphi_i^{(m)} \,dx = \int_{\Omega} \left( \frac{\eta_p}{\lambda} (c^{(m)} - I_d) \right) : \nabla^{\mathcal{M}} \varphi^{(m)} \,dx,$$

and the results follows thanks to the convergence of $\nabla^{\mathcal{M}}\varphi^{(m)}(x)$ to $\nabla\varphi(x)$ for almost every $x \in \Omega$. $\square$

This result extends to nonlinear laws, under stronger convergence and boundedness assumptions (and/or regularization of the function $g$) for the sequence of conformation tensors (in fine, what is needed is the weak convergence in $L^1(\Omega)^{d \times d}$ of $g(c^{(m)})(c^{(m)} - I_d)$).
An easy calculation shows that this relation may be recast as a finite volume formulation, in the sense that the right-hand side may be seen as a sum over the faces of $K_{i-\frac{1}{2},j}^x$ of a discretization of the flux associated with $\text{div}\tau_p$, i.e. the integral of the first component of $\tau_p n_{K,\sigma}$ (see Fig. 3).

- On the vertical face $\varepsilon = \varepsilon_{i,j}^x$, the flux is given by

$$ \int_{\varepsilon} (\tau_p n_{\varepsilon})^x = \frac{1}{h_j^y} (\tau_p)_{i,j}^{xx} = |\varepsilon_{i,j}^x| (\tau_p)_{i,j}^{xx}. $$

- The flux through the horizontal half face $\varepsilon = \varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x \cap \bar{K}_{i,j}$ reads

$$ \int_{\varepsilon} (\tau_p n_{\varepsilon})^x = \frac{1}{h_{j+\frac{1}{2}}^y} \left( \frac{|K_{i,j}|}{4} (\tau_p)_{i,j}^{xy} + \frac{|K_{i,j+1}|}{4} (\tau_p)_{i,j+1}^{xy} \right) = |\varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x| \left( \alpha (\tau_p)_{i,j}^{xy} + (1-\alpha) (\tau_p)_{i,j+1}^{xy} \right), $$

with

$$ \alpha = \frac{h_j^y}{2h_{j+\frac{1}{2}}^y}, \quad \text{and} \quad 1-\alpha = \frac{h_{j+1}^y}{2h_{j+\frac{1}{2}}^y}. $$

Hence, as usual when such a duality technique is used, the approximation of the tensor at the horizontal faces may seem strange. It is a convex combination of the unknowns in the two neighbouring cells, but with coefficients that are not those which would be given by a linear interpolation. Indeed, on $\varepsilon = \varepsilon_{i-\frac{1}{4},j+\frac{1}{2}}^x$, we have

$$ (\tau_p)_\varepsilon^{xy} = \alpha (\tau_p)_{i,j}^{xy} + (1-\alpha) (\tau_p)_{i,j+1}^{xy}, $$

whereas, with a linear interpolation, we would have

$$ (\tau_p)_\varepsilon^{xy} = (1-\alpha) (\tau_p)_{i,j}^{xy} + \alpha (\tau_p)_{i,j+1}^{xy}, $$

so that the formulae differ on a non-uniform mesh.

Up to now, we have only invoked stability arguments to justify the derivation of the discrete divergence of the elastic extra-stress tensor through a weak formulation. However, thanks to the consistency properties of the discrete gradient $\nabla^M$, this technique also ensures by construction the weak consistency of this discretization in the Lax-Wendroff sense, as stated in the following result for the Oldroyd-B model. For short, it is written for functions depending only on space but trivially extends to time-dependent functions.

**Lemma 3.2.** Let $(\mathcal{M}^{(m)})_{m \in \mathbb{N}}$ be sequence of structured discretizations of $\Omega$, with a space step tending to zero. Let $(c^{(m)})_{m \in \mathbb{N}}$ be a sequence of discrete conformation tensors, with $c^{(m)}$ associated with $\mathcal{M}^{(m)}$. We suppose that $(c^{(m)})_{m \in \mathbb{N}}$ weakly converges in $L^1(\Omega)^{d \times d}$ to $c$. Let $\varphi \in C_c^\infty(\Omega)^d$, and let us denote by $\varphi^{(m)}$ the interpolate of $\varphi$ in the space of velocity discrete functions defined by $\varphi^{(m)}(x) = \varphi_i(x_\sigma)$ for $x \in D_\sigma$, $\sigma \in (\varepsilon^{(m)})^{(i)}$, $1 \le i \le d$, with $x_\sigma$ the mass center of $\sigma$. Finally, for a discrete tensor $T$, i.e. a tensor-valued piecewise constant function over the primal cells of $\mathcal{M}^{(m)}$, we denote by $\text{div}_{\mathcal{M}^{(m)}} T: \Omega \to \mathbb{R}^d$ the vector-valued function defined by $(\text{div}_{\mathcal{M}^{(m)}} T)_i = \text{div}_{\sigma,i} T$ for $x \in D_\sigma$, $\sigma \in (\varepsilon^{(m)})^{(i)}$, $1 \le i \le d$. Then

$$ \lim_{m \to +\infty} \int_{\Omega} \text{div}_{\mathcal{M}^{(m)}} \left( \frac{\eta_p}{\lambda} (c^{(m)} - I_d) \right) \cdot \varphi_i^{(m)} \, dx = - \int_{\Omega} \frac{\eta_p}{\lambda} (c - I_d) : \nabla \varphi \, dx. $$

*Proof.* By construction, we have

$$ \int_{\Omega} \text{div}_{\mathcal{M}^{(m)}} \left( \frac{\eta_p}{\lambda} (\mathbf{c}^{(m)} - \mathbf{I}_d) \right) \cdot \varphi_i^{(m)} \, d\mathbf{x} = \int_{\Omega} \left( \frac{\eta_p}{\lambda} (\mathbf{c}^{(m)} - \mathbf{I}_d) \right) : \nabla^{\mathcal{M}} \varphi^{(m)} \, d\mathbf{x}, $$

and the results follows thanks to the convergence of $\nabla^{\mathcal{M}}\varphi^{(m)}(x)$ to $\nabla\varphi(x)$ for almost every $x \in \Omega$. $\square$

This result extends to nonlinear laws, under stronger convergence and boundedness assumptions (and/or regularization of the function $g$) for the sequence of conformation tensors (in fine, what is needed is the weak convergence in $L^1(\Omega)^{d \times d}$ of $g(c^{(m)})(c^{(m)} - \mathbf{I}_d)$).
## 4. A MORE EFFICIENT AND ACCURATE PRESSURE CORRECTION SCHEME

Unfortunately, computations show that the scheme (21) is not as accurate and efficient as expected. The reasons of this behaviour may be traced back and cured as follows.

- First, the update of the conformation tensor in the projection step (i.e. the term $\text{div}_{\sigma,i} \tau_p^{n+1} - \text{div}_{\sigma,i} \tau_p^n$) tightly couples the equations of this step and thus makes its resolution time-consuming. In addition, it does not seem to bring any significant gain in accuracy. This term is thus suppressed, so as to decouple the transport equation of the conformation tensor (21d) from equations (21b) and (21c) (which, as usual, are recast under the form of an elliptic problem for the pressure).

- Second, especially at high Weissenberg numbers, the production terms in equation (21d) (i.e. the terms $\nabla u_K^{n+1} c_K^{n+1} + (c_K^{n+1})^t (\nabla u_K^{n+1})^t$) are likely to become very stiff, so that stability and accuracy of usual schemes may require prohibitively small space and time steps. The analysis of a toy scalar problem (see Sect. D) shows that the scheme accuracy is dramatically enhanced by making a change of variable from $c$ to $\log(c)$. However, for general equations of state (in particular, as in the FENE-CR model, when the function $g(c)$ does not boil down to identity), this change of variable brings additional nonlinearities. A solution for this problem is given by the observation that, still on the toy problem, the scheme accuracy is essentially preserved by a 3-steps Strang-type decoupling, that consists in performing one half-step of homogeneous transport of $\log(c)$, then dealing with the reaction terms with the original variable $c$ and finishing by the second half-step of transport of $\log(c)$. Up to the choice of the variable in the middle step, a similar strategy was already adopted in [43], with a different space discretization. With the specific space discretization considered here, i.e. with piecewise-constant discrete conformation tensors, the middle step consists in a family of Ordinary Differential equations, one per primal cell, for all the tensor components. Since the stiffness of this step is directly correlated to the magnitude of the velocity gradient in the cell, this allows us to adopt a local time-stepping strategy with a time step per cell, which dramatically enhances the accuracy and robustness of the scheme. Moreover, in most applications, we observe that large velocity gradients are localized in very specific zones of the flow so that small local time steps are used only in a few cells and the scheme efficiency is only slightly affected.

Finally, we obtain a two-steps scheme, where we successively deal first with hydrodynamics and then with the conformation tensor transport. The hydrodynamics step reads

**Prediction step** – Solve for $\tilde{u}^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \\ \qquad \frac{\rho}{\delta t} (\tilde{u}_{\sigma,i}^{n+1} - u_{\sigma,i}^n) + \xi \rho \text{div}_\sigma (\tilde{u}_i^{n+1} u^n) - \text{div}_{\sigma,i} \tau_s (\tilde{u}^{n+1}) + \nabla_{\sigma,i} p^n - \text{div}_{\sigma,i} \tau_p^n = 0. \quad (32a) $$

**Correction step** – Solve for $p^{n+1}$ and $u^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \quad \frac{\rho}{\delta t} (u_{\sigma,i}^{n+1} - \tilde{u}_{\sigma,i}^{n+1}) + \nabla_{\sigma,i} (p^{n+1} - p^n) = 0, \quad (32b) $$

$$ \forall K \in \mathcal{M}, \quad \text{div}_K (u^{n+1}) = 0. \quad (32c) $$

Up to the (explicit) term $\text{div}_{\sigma,i} \tau_p^n$ in the prediction step, this system is the usual so-called incremental pressure-correction scheme for incompressible flows. All the discrete operators featured in these equations have already been described in Section 3.

The discretization of the constitutive equation (1c) is then performed with the following “Strang-log” scheme

**Advection I** – Solve for $c^{n+\frac{1}{3}}$:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t/2} \left( \log c_K^{n+\frac{1}{3}} - \log c_K^n \right) + \text{div}_K \left( u^{n+1} \log c^{n+\frac{1}{3}} \right) = 0. \quad (33a) $$
## 4. A MORE EFFICIENT AND ACCURATE PRESSURE CORRECTION SCHEME

Unfortunately, computations show that the scheme (21) is not as accurate and efficient as expected. The reasons of this behaviour may be traced back and cured as follows.

- First, the update of the conformation tensor in the projection step (i.e. the term $\operatorname{div}_{\sigma,i} \boldsymbol{\tau}_p^{n+1} - \operatorname{div}_{\sigma,i} \boldsymbol{\tau}_p^n$) tightly couples the equations of this step and thus makes its resolution time-consuming. In addition, it does not seem to bring any significant gain in accuracy. This term is thus suppressed, so as to decouple the transport equation of the conformation tensor (21d) from equations (21b) and (21c) (which, as usual, are recast under the form of an elliptic problem for the pressure).

- Second, especially at high Weissenberg numbers, the production terms in equation (21d) (i.e. the terms $\nabla u_K^{n+1} c_K^{n+1} + (c_K^{n+1})^t (\nabla u_K^{n+1})^t$) are likely to become very stiff, so that stability and accuracy of usual schemes may require prohibitively small space and time steps. The analysis of a toy scalar problem (see Sect. D) shows that the scheme accuracy is dramatically enhanced by making a change of variable from $c$ to $\log(c)$. However, for general equations of state (in particular, as in the FENE-CR model, when the function $g(c)$ does not boil down to identity), this change of variable brings additional nonlinearities. A solution for this problem is given by the observation that, still on the toy problem, the scheme accuracy is essentially preserved by a 3-steps Strang-type decoupling, that consists in performing one half-step of homogeneous transport of $\log(c)$, then dealing with the reaction terms with the original variable $c$ and finishing by the second half-step of transport of $\log(c)$. Up to the choice of the variable in the middle step, a similar strategy was already adopted in [43], with a different space discretization. With the specific space discretization considered here, i.e. with piecewise-constant discrete conformation tensors, the middle step consists in a family of Ordinary Differential equations, one per primal cell, for all the tensor components. Since the stiffness of this step is directly correlated to the magnitude of the velocity gradient in the cell, this allows us to adopt a local time-stepping strategy with a time step per cell, which dramatically enhances the accuracy and robustness of the scheme. Moreover, in most applications, we observe that large velocity gradients are localized in very specific zones of the flow so that small local time steps are used only in a few cells and the scheme efficiency is only slightly affected.

Finally, we obtain a two-steps scheme, where we successively deal first with hydrodynamics and then with the conformation tensor transport. The hydrodynamics step reads

**Prediction step** – Solve for $\tilde{u}^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \\ \qquad \frac{\rho}{\delta t} (\tilde{u}_{\sigma,i}^{n+1} - u_{\sigma,i}^n) + \xi \rho \operatorname{div}_\sigma (\tilde{u}_i^{n+1} u^n) - \operatorname{div}_{\sigma,i} \boldsymbol{\tau}_s (\tilde{u}^{n+1}) + \nabla_{\sigma,i} p^n - \operatorname{div}_{\sigma,i} \boldsymbol{\tau}_p^n = 0. \quad (32a) $$

**Correction step** – Solve for $p^{n+1}$ and $u^{n+1}$:

$$ \text{For } 1 \le i \le d, \quad \forall \sigma \in \mathcal{E}_S^{(i)}, \quad \frac{\rho}{\delta t} (u_{\sigma,i}^{n+1} - \tilde{u}_{\sigma,i}^{n+1}) + \nabla_{\sigma,i} (p^{n+1} - p^n) = 0, \quad (32b) $$

$$ \forall K \in \mathcal{M}, \quad \operatorname{div}_K (u^{n+1}) = 0. \quad (32c) $$

Up to the (explicit) term $\operatorname{div}_{\sigma,i} \boldsymbol{\tau}_p^n$ in the prediction step, this system is the usual so-called incremental pressure-correction scheme for incompressible flows. All the discrete operators featured in these equations have already been described in Section 3.

The discretization of the constitutive equation (1c) is then performed with the following “Strang-log” scheme

**Advection I** – Solve for $c^{n+\frac{1}{3}}$:

$$ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t/2} \left( \log c_K^{n+\frac{1}{3}} - \log c_K^n \right) + \operatorname{div}_K \left( u^{n+1} \log c^{n+\frac{1}{3}} \right) = 0. \quad (33a) $$
$$ \begin{gathered} \text{ODE - Set } \mathbf{c}_K(t_n) = \mathbf{c}_K^{n+\frac{1}{3}} \text{ and solve for } \mathbf{c}^{n+\frac{2}{3}} = \mathbf{c}_K(t_n + \delta t): \\ \forall K \in \mathcal{M}, \quad \partial_t \mathbf{c}_K - (\nabla_K \mathbf{u}^{n+1}) \mathbf{c}_K - \mathbf{c}_K^t (\nabla_K \mathbf{u}^{n+1})^t + \frac{g(\mathbf{c}_K)}{\lambda} (\mathbf{c}_K - \mathbf{I}_d) = 0. \quad (33b) \end{gathered} $$

$$ \begin{gathered} \text{Advection II - Solve for } \mathbf{c}^{n+1}: \\ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t/2} \left( \log \mathbf{c}_K^{n+1} - \log \mathbf{c}_K^{n+\frac{2}{3}} \right) + \text{div}_K (\mathbf{u}^{n+1} \log \mathbf{c}^{n+1}) = 0. \quad (33c) \end{gathered} $$

Transport steps are discretized by a standard backward first-order upwind scheme. From the form of these equations, these steps preserve the symmetry of the conformation tensor: if $\mathbf{c}^n$ is symmetric, the same equation is solved for $\mathbf{c}_{i,j}^{n+1}$ and $\mathbf{c}_{j,i}^{n+1}$, $1 \le i, j \le d$. In addition, since we search for $\log(\mathbf{c})$ and then compute $\mathbf{c}$ by taking the exponential of the result, these transport steps necessarily yields positive definite conformation tensors. The local ODE (33b) is solved using a first-order semi-implicit Euler scheme, with a local sub-time step $\delta t_K = \delta t/N_K$, that reads, for $K \in \mathcal{M}$ and $0 \le k \le N_K$

$$ \frac{1}{\delta t_K} (\mathbf{c}_K^{k+1} - \mathbf{c}_K^k) = \nabla \mathbf{u}_K^{n+1} \mathbf{c}_K^{k+1} + (\mathbf{c}_K^{k+1})^t (\nabla \mathbf{u}_K^{n+1})^t - \frac{g(\mathbf{c}_K^k)}{\lambda} (\mathbf{c}_K^{k+1} - \mathbf{I}_d), \quad (34) $$

with $\mathbf{c}_K^0 = \mathbf{c}_K^{n+\frac{1}{3}}$ and $\mathbf{c}_K^{n+\frac{2}{3}} = \mathbf{c}_K^{N_K}$. We prove in Section A.2 that this relation yields a symmetric positive definite tensor under the constraint for the choice of the sub-time step $\delta t_K$ stated by Inequality (A.5) of Appendix A.2, which imposes a scaling of the time step as $1/\|\nabla \mathbf{u}_K^{n+1}\|_{\infty}$. In addition, for the FENE-CR model, we prove in Section B that, under additional constraints on the sub-time step (see Lem. B.2), the solution of (34) satisfies $\text{tr}(\mathbf{c}_K^{k+1}) < b$, for $K \in \mathcal{M}$.

Note that the constitutive relation only appears in the step (33b) and in the expression of the polymeric stress tensor in the prediction step (32a), with explicit-in-time values in the latter case; in addition, it associates a piecewise constant conformation tensor to a piecewise constant (on the same mesh) stress tensor. These two features make changes of constitutive relations very easy in practice. The ODE step (34) is already written under its general form and, for the velocity prediction step (32a), the tensor $\boldsymbol{\tau}_p^n$ reads:

$$ (\boldsymbol{\tau}_p^n)_K = \frac{\eta_p}{\lambda} \mathbf{f}(\mathbf{c}_K^n) (\mathbf{c}_K^n - \mathbf{I}_d), \quad \text{for } K \in \mathcal{M}, $$

and the discretization of its divergence is described in Section 3.3.1 (RT discretization) or Section 3.3.2 (MAC discretization).

## 5. NUMERICAL RESULTS

Here we consider the dimensionless version of the system of equation (1) in the vanishing inertia limit ($\xi = 0$),

$$ \text{Re} \partial_t \mathbf{u} = -\nabla p + \beta \text{div} (\nabla \mathbf{u} + (\nabla \mathbf{u})^t) + \frac{1-\beta}{\text{W}_i} \text{div} (\mathbf{f}(\mathbf{c})(\mathbf{c}-\mathbf{I}_d)), \quad (35a) $$

$$ \text{div } \mathbf{u} = 0, \quad (35b) $$

$$ \partial_t \mathbf{c} + \mathbf{u} \cdot \nabla \mathbf{c} - (\nabla \mathbf{u}) \mathbf{c} - \mathbf{c} (\nabla \mathbf{u})^t + \frac{1}{\text{W}_i} \mathbf{f}(\mathbf{c})(\mathbf{c} - \mathbf{I}_d) = 0, \quad (35c) $$

where space is non-dimensionalized with $L$ and time with $L/U$, where $L$ and $U$ stand for the characteristic length scale and velocity respectively, $\text{W}_i = \lambda U/L$, $\beta = \eta_s/(\eta_s + \eta_p)$ is the viscosity ratio and Re is the Reynolds number defined by $\text{Re} = \rho U L/(\eta_s + \eta_p)$. The pressure is non-dimensionalized with a reference pressure equal to $(\eta_s + \eta_p) U/L$. The constitutive relation is given by

$$ \mathbf{f}(\mathbf{c}) = \mathbf{I}_d \text{ for the Oldroyd-B model, } \mathbf{f}(\mathbf{c}) = \frac{b}{b - \text{tr}(\mathbf{c})} \mathbf{I}_d, b > d, \text{ for FENE-CR.} $$
$$ \begin{gathered} \text{ODE - Set } \mathbf{c}_K(t_n) = \mathbf{c}_K^{n+\frac{1}{3}} \text{ and solve for } \mathbf{c}^{n+\frac{2}{3}} = \mathbf{c}_K(t_n + \delta t): \\ \forall K \in \mathcal{M}, \quad \partial_t \mathbf{c}_K - (\nabla_K \mathbf{u}^{n+1}) \mathbf{c}_K - \mathbf{c}_K^t (\nabla_K \mathbf{u}^{n+1})^t + \frac{g(\mathbf{c}_K)}{\lambda} (\mathbf{c}_K - \mathbf{I}_d) = 0. \quad (33b) \end{gathered} $$

$$ \begin{gathered} \text{Advection II - Solve for } \mathbf{c}^{n+1}: \\ \forall K \in \mathcal{M}, \quad \frac{1}{\delta t/2} \left( \log \mathbf{c}_K^{n+1} - \log \mathbf{c}_K^{n+\frac{2}{3}} \right) + \text{div}_K (\mathbf{u}^{n+1} \log \mathbf{c}^{n+1}) = 0. \quad (33c) \end{gathered} $$

Transport steps are discretized by a standard backward first-order upwind scheme. From the form of these equations, these steps preserve the symmetry of the conformation tensor: if $\mathbf{c}^n$ is symmetric, the same equation is solved for $\mathbf{c}_{i,j}^{n+1}$ and $\mathbf{c}_{j,i}^{n+1}$, $1 \le i, j \le d$. In addition, since we search for $\log(\mathbf{c})$ and then compute $\mathbf{c}$ by taking the exponential of the result, these transport steps necessarily yields positive definite conformation tensors. The local ODE (33b) is solved using a first-order semi-implicit Euler scheme, with a local sub-time step $\delta t_K = \delta t/N_K$, that reads, for $K \in \mathcal{M}$ and $0 \le k \le N_K$

$$ \frac{1}{\delta t_K} (\mathbf{c}_K^{k+1} - \mathbf{c}_K^k) = \nabla \mathbf{u}_K^{n+1} \mathbf{c}_K^{k+1} + (\mathbf{c}_K^{k+1})^t (\nabla \mathbf{u}_K^{n+1})^t - \frac{g(\mathbf{c}_K^k)}{\lambda} (\mathbf{c}_K^{k+1} - \mathbf{I}_d), \quad (34) $$

with $\mathbf{c}_K^0 = \mathbf{c}_K^{n+\frac{1}{3}}$ and $\mathbf{c}_K^{n+\frac{2}{3}} = \mathbf{c}_K^{N_K}$. We prove in Section A.2 that this relation yields a symmetric positive definite tensor under the constraint for the choice of the sub-time step $\delta t_K$ stated by Inequality (A.5) of Appendix A.2, which imposes a scaling of the time step as $1/\|\nabla \mathbf{u}_K^{n+1}\|_{\infty}$. In addition, for the FENE-CR model, we prove in Section B that, under additional constraints on the sub-time step (see Lem. B.2), the solution of (34) satisfies $\text{tr}(\mathbf{c}_K^{k+1}) < b$, for $K \in \mathcal{M}$.

Note that the constitutive relation only appears in the step (33b) and in the expression of the polymeric stress tensor in the prediction step (32a), with explicit-in-time values in the latter case; in addition, it associates a piecewise constant conformation tensor to a piecewise constant (on the same mesh) stress tensor. These two features make changes of constitutive relations very easy in practice. The ODE step (34) is already written under its general form and, for the velocity prediction step (32a), the tensor $\boldsymbol{\tau}_p^n$ reads:

$$ (\boldsymbol{\tau}_p^n)_K = \frac{\eta_p}{\lambda} \mathbf{f}(\mathbf{c}_K^n) (\mathbf{c}_K^n - \mathbf{I}_d), \quad \text{for } K \in \mathcal{M}, $$

and the discretization of its divergence is described in Section 3.3.1 (RT discretization) or Section 3.3.2 (MAC discretization).

## 5. NUMERICAL RESULTS

Here we consider the dimensionless version of the system of equation (1) in the vanishing inertia limit ($\xi = 0$),

$$ \text{Re} \partial_t \mathbf{u} = -\nabla p + \beta \text{div} (\nabla \mathbf{u} + (\nabla \mathbf{u})^t) + \frac{1-\beta}{\text{W}_i} \text{div} (\mathbf{f}(\mathbf{c})(\mathbf{c}-\mathbf{I}_d)), \quad (35a) $$

$$ \text{div} \mathbf{u} = 0, \quad (35b) $$

$$ \partial_t \mathbf{c} + \mathbf{u} \cdot \nabla \mathbf{c} - (\nabla \mathbf{u}) \mathbf{c} - \mathbf{c} (\nabla \mathbf{u})^t + \frac{1}{\text{W}_i} \mathbf{f}(\mathbf{c})(\mathbf{c} - \mathbf{I}_d) = 0, \quad (35c) $$

where space is non-dimensionalized with $L$ and time with $L/U$, where $L$ and $U$ stand for the characteristic length scale and velocity respectively, $\text{W}_i = \lambda U/L$, $\beta = \eta_s/(\eta_s + \eta_p)$ is the viscosity ratio and Re is the Reynolds number defined by $\text{Re} = \rho UL/(\eta_s + \eta_p)$. The pressure is non-dimensionalized with a reference pressure equal to $(\eta_s + \eta_p) U/L$. The constitutive relation is given by

$$ \mathbf{f}(\mathbf{c}) = \mathbf{I}_d \text{ for the Oldroyd-B model, } \mathbf{f}(\mathbf{c}) = \frac{b}{b - \text{tr}(\mathbf{c})} \mathbf{I}_d, b > d, \text{ for FENE-CR.} $$
The momentum balance equation (35a) may be complemented by a body force; in this case, it is explicitly mentioned. In all this section, we solve equation (35) with $Re = 1$ (remember that $Re$ multiplies the velocity time derivative term only, the momentum convection term being removed in (35)), so the data of $\beta$, $Wi$, the geometry and the initial and boundary conditions completely defines the problem.

We address two benchmark cases from the literature. The first one is the lid-driven cavity, which we solve with the MAC scheme and with the RT discretization, both on Cartesian (uniform and non-uniform) grids. The second one is the flow around a cylinder confined in a channel, which we solve on an unstructured body-fitting mesh using the RT discretization. Even though classic [18, 20, 22, 23, 27, 32, 43, 48], each of these cases is challenging to solve numerically even at moderate $Wi$ numbers. In the case of the lid-driven cavity, it is only with the recent introduction of the log-conformation formulation that results at $Wi \ge 1$ have been obtained [22, 23]. For the flow around the cylinder, a singularity in the wake of the cylinder appears at $Wi \approx 0.7$ that is very difficult to capture accurately [20, 32].

Our goals are here to validate the proposed scheme by comparison with results from the literature, and to study its accuracy and numerical convergence properties. We present results for both the Oldroyd-B and the FENE-CR models, these results being new in the latter case, and compare the obtained solutions. In the lid-driven cavity case, we quantitatively assess the convergence of the solutions of the FENE-CR model to the solutions of the Oldroyd-B model when the parameter $b$ limiting the trace of the conformation tensor tends to $+\infty$. To perform these computations, the proposed scheme was implemented in the open-source CALIF³S software developed by the Institut de Radioprotection et de Sûreté Nucléaire IRSN [16]. The whole domain is decomposed into several sub-domains with the METIS (4.0.3) graph partitioner and the communication between the sub-domains is done through the OpenMPI (3.1.5) interface. The linear systems are solved by Krylov methods, using a Block-Jacobi preconditioner at the global level and an ILU preconditioner for the blocks associated to the subdomains. Parallel computations were then carried out on TotalEnergies' group supercomputer Pangea II.

## 5.1. Lid-driven cavity

In this test, the computational domain is the unit-square $\Omega = (0, 1)^2$. At the initial time, the fluid is at rest (i.e. $u(x, 0) = 0$ a.e. in $\Omega$) and the conformation tensor is uniformly equal to the identity ($c(x, 0) = I_d$ a.e. in $\Omega$). Dirichlet boundary conditions are imposed to the velocity components at each side of the square, with $u = (8x^2(1-x)^2(1+\tanh(8t-4)), 0)$ (see [22]) on the lid and $u = 0$ on the three other sides and at any time. The viscosity ratio is $\beta = 0.5$.

We solve this problem for a sequence of successively refined Cartesian meshes. The three coarsest meshes are uniform $64 \times 64$, $128 \times 128$ and $256 \times 256$ grids. The other ones, denoted by $Mn, n = 1, \dots, 4$, use a non-uniform step in the two directions. These are built from lists for the $x$-coordinates and $y$-coordinates of the edges, which are defined as follows.

- The set of $x$-coordinates is $\mathcal{G}((0, 0.001/n), (0.5, 0.004/n)) \cup \mathcal{G}((0.5, 0.004/n), (1, 0.001/n))$,

- The set of $y$-coordinates is $\mathcal{G}((0, 0.004/n), (1, 0.001/n))$, where $\mathcal{G}((a, \delta x_i), (b, \delta x_f))$ is the set of real numbers beginning at $a$, ending at $b$, with a step following a geometrical progression from $\delta x_i$ at $a$ to $\delta x_f$ at $b$. With this procedure, the number of cells for the M4 mesh is close to 5.2 millions.

In these computations, even if, as we will see, a steady state solution is obtained, the time step must be monitored carefully to ensure time convergence. If the time step is too large, the solution shows a spurious unstationary behaviour. Hence, unless explicitly mentioned, the time step is $\delta t = 10^{-5}$ at the initial time and then increases, but is limited to keep the relative difference between the beginning and end-of step values of the velocity and the conformation tensor components lower than 0.001 in the maximum norm. The time step is also not allowed to be greater than $\delta t = 10^{-3}$, which is the time step beyond which unstationary phenomena appear. The sub-time-step for the solution of the ODE on the cell $K$ is set to $\delta t/n_e$ with $n_e$ the smallest integer number such that $\delta t/n_e \le 1/(2m ||\nabla_K u^{n+1}||_\infty)$ for the Oldroyd-B model and $\delta t/n_e \le (b-2)/(4\lambda bm ||\nabla_K u^{n+1}||_\infty^2)$ for the FENE-CR model, with $m = 100$. These bounds are here to preserve the positive definiteness of the
The momentum balance equation (35a) may be complemented by a body force; in this case, it is explicitly mentioned. In all this section, we solve equation (35) with $Re = 1$ (remember that $Re$ multiplies the velocity time derivative term only, the momentum convection term being removed in (35)), so the data of $\beta$, $Wi$, the geometry and the initial and boundary conditions completely defines the problem.

We address two benchmark cases from the literature. The first one is the lid-driven cavity, which we solve with the MAC scheme and with the RT discretization, both on Cartesian (uniform and non-uniform) grids. The second one is the flow around a cylinder confined in a channel, which we solve on an unstructured body-fitting mesh using the RT discretization. Even though classic [18, 20, 22, 23, 27, 32, 43, 48], each of these cases is challenging to solve numerically even at moderate $Wi$ numbers. In the case of the lid-driven cavity, it is only with the recent introduction of the log-conformation formulation that results at $Wi \ge 1$ have been obtained [22, 23]. For the flow around the cylinder, a singularity in the wake of the cylinder appears at $Wi \approx 0.7$ that is very difficult to capture accurately [20, 32].

Our goals are here to validate the proposed scheme by comparison with results from the literature, and to study its accuracy and numerical convergence properties. We present results for both the Oldroyd-B and the FENE-CR models, these results being new in the latter case, and compare the obtained solutions. In the lid-driven cavity case, we quantitatively assess the convergence of the solutions of the FENE-CR model to the solutions of the Oldroyd-B model when the parameter $b$ limiting the trace of the conformation tensor tends to $+\infty$. To perform these computations, the proposed scheme was implemented in the open-source CALIF³S software developed by the Institut de Radioprotection et de Sûreté Nucléaire IRSN [16]. The whole domain is decomposed into several sub-domains with the METIS (4.0.3) graph partitioner and the communication between the sub-domains is done through the OpenMPI (3.1.5) interface. The linear systems are solved by Krylov methods, using a Block-Jacobi preconditioner at the global level and an ILU preconditioner for the blocks associated to the subdomains. Parallel computations were then carried out on TotalEnergies' group supercomputer Pangea II.

## 5.1. Lid-driven cavity

In this test, the computational domain is the unit-square $\Omega = (0, 1)^2$. At the initial time, the fluid is at rest (i.e. $u(x, 0) = 0$ a.e. in $\Omega$) and the conformation tensor is uniformly equal to the identity ($c(x, 0) = I_d$ a.e. in $\Omega$). Dirichlet boundary conditions are imposed to the velocity components at each side of the square, with $u = (8x^2(1-x)^2(1+\tanh(8t-4)),$ 0) (see [22]) on the lid and $u = 0$ on the three other sides and at any time. The viscosity ratio is $\beta = 0.5$.

We solve this problem for a sequence of successively refined Cartesian meshes. The three coarsest meshes are uniform $64 \times 64$, $128 \times 128$ and $256 \times 256$ grids. The other ones, denoted by $Mn, n = 1, \dots, 4$, use a non-uniform step in the two directions. These are built from lists for the $x$-coordinates and $y$-coordinates of the edges, which are defined as follows.

- The set of $x$-coordinates is $\mathcal{G}((0, 0.001/n), (0.5, 0.004/n)) \cup \mathcal{G}((0.5, 0.004/n), (1, 0.001/n))$,

- The set of $y$-coordinates is $\mathcal{G}((0, 0.004/n), (1, 0.001/n))$, where $\mathcal{G}((a, \delta x_i), (b, \delta x_f))$ is the set of real numbers beginning at $a$, ending at $b$, with a step following a geometrical progression from $\delta x_i$ at $a$ to $\delta x_f$ at $b$. With this procedure, the number of cells for the M4 mesh is close to 5.2 millions.

In these computations, even if, as we will see, a steady state solution is obtained, the time step must be monitored carefully to ensure time convergence. If the time step is too large, the solution shows a spurious unstationary behaviour. Hence, unless explicitly mentioned, the time step is $\delta t = 10^{-5}$ at the initial time and then increases, but is limited to keep the relative difference between the beginning and end-of step values of the velocity and the conformation tensor components lower than 0.001 in the maximum norm. The time step is also not allowed to be greater than $\delta t = 10^{-3}$, which is the time step beyond which unstationary phenomena appear. The sub-time-step for the solution of the ODE on the cell $K$ is set to $\delta t/n_e$ with $n_e$ the smallest integer number such that $\delta t/n_e \le 1/(2m ||\nabla_K u^{n+1}||_\infty)$ for the Oldroyd-B model and $\delta t/n_e \le (b-2)/(4\lambda bm ||\nabla_K u^{n+1}||_\infty^2)$ for the FENE-CR model, with $m = 100$. These bounds are here to preserve the positive definiteness of the
conformation tensor when solving the ODE by a backward Euler scheme, see Sections A.2 and B. The CPU-time used in this step remains almost negligible (less than 3% of the total time), so a more sophisticated algorithm would not significantly increase the efficiency.

### 5.1.1. Oldroyd-B model at $Wi = 1$ and $Wi = 3$

For these two test cases, we observe the same convergence behaviour for the MAC and RT discretizations. For the finest meshes, results cannot be visually distinguished, except for the specific case of the conformation tensor at the cavity lid (see the discussion below). We present here the results obtained with the MAC scheme.

We start by plotting in Figure 6 the 2D representations of the velocity magnitude and the trace of the conformation tensor for $Wi = 1$ and $Wi = 3$ when the steady state is reached. The velocity field exhibits some of the classical features that have been observed in the literature. Each corner at the bottom of the cavity contains a vortex at very low velocity that effectively behaves as a stagnant zone. As we increase the $Wi$, we transition from a left-right symmetry ($Wi = 0$) to an asymmetric profile with a shift of the primary central vortex and of the region of maximum velocity close to the lid towards the left-hand side. The most interesting features are observed for the trace of the conformations tensor. We see a dramatic increase of the maximum value of the trace from $Wi = 1$ to $Wi = 3$, going from approximately $5 \times 10^3$ to more than $10^5$. This maximum is reached at the lid with zones of large stresses both in the middle of the lid and in the top right corner. This suggests that mesh convergence in these regions may be particularly difficult to obtain. We also observe that regions of high stresses have extended further inside the cavity at $Wi = 3$ compared to $Wi = 1$, which may impact convergence even away from the lid.

We then study mesh convergence and compare quantitatively our scheme to results from the literature, starting with the Weissenberg number equal to 1 and comparing to results in [23]. The computation is performed up to a time $t = 10$, with a steady state reached for $t \approx 8$. The first component of the velocity, $u_x$, along the line $x = 0.5$, is presented in Figure 7a. The steady state is almost independent from the mesh and is in close agreement with results in [23]. The second component of the velocity, $u_y$, along the line $y = 0.75$ is plotted in Figure 7b. Convergence is slower and seems to be roughly obtained for the $256 \times 256$ grid, which is also in close agreement with [23].

The most difficult aspect of these computations is to accurately estimate the conformation tensor near the lid. We investigate this issue here using the $Mn$ meshes that are specifically designed for this purpose, with finer meshes in the region close to the lid. Figure 7c shows $\text{tr}(c)$ along the line $x = 0.5$ for M1, M2, M3 and M4. On this plot, convergence seems to be achieved for any of the grids. However, we can see in Figure 7d that the profile of $\text{tr}(c)$ becomes very steep at the top. Figures 7e and 7f show $\text{tr}(c)$ along the lines $y = 0.975$ and $y = 1$ respectively, for the different meshes. At $y = 0.975$, convergence seems almost achieved. The picture is completely different at $y = 1$ where the maximum value of $\text{tr}(c)$, obtained close to $x = 0.5$, increases dramatically when refining the mesh. Convergence is very slow as, even for the finest mesh, we see a multiplication of the value of the maximum by a factor of about 1.5 each time the mesh step is divided by $(n+1)/n$ for $n = 1, \dots, 3$.

We now focus on the case $Wi = 3$ and compare with results in [18,23]. When the mesh is coarse, an unsteady behavior is observed that disappears when refining the mesh. The temporal evolution of the kinetic energy,

$$E_{\text{kin}}(t) = \frac{1}{2} \int_{\Omega} |\mathbf{u}(\mathbf{x}, t)|^2 d\mathbf{x},$$

is presented in Figure 8. For the uniform $64 \times 64$ grid, we observe periodic oscillations of the kinetic energy that are related to the appearance of vortices in the upper-right corner of the cavity. However, this phenomenon seems to be purely numerical, since, with finer meshes, oscillations disappear and a stationary solution is obtained after $t \approx 40$. Such numerical instabilities had also been observed in [23], although they did not obtain a stable solution by refining the mesh. In [18], solutions for this case may also be slightly unsteady. We compare our stationary solution to the solutions obtained in these two papers at a given time $t = 40$. Figures 9a and 9b depict profiles of the velocity field for the $Mn$ meshes, showing that convergence is slower than for $Wi = 1$. We find that our results for the first velocity component (Fig. 9a) are in very good agreement with those of
conformation tensor when solving the ODE by a backward Euler scheme, see Sections A.2 and B. The CPU-time used in this step remains almost negligible (less than 3% of the total time), so a more sophisticated algorithm would not significantly increase the efficiency.

### 5.1.1. Oldroyd-B model at $Wi = 1$ and $Wi = 3$

For these two test cases, we observe the same convergence behaviour for the MAC and RT discretizations. For the finest meshes, results cannot be visually distinguished, except for the specific case of the conformation tensor at the cavity lid (see the discussion below). We present here the results obtained with the MAC scheme.

We start by plotting in Figure 6 the 2D representations of the velocity magnitude and the trace of the conformation tensor for $Wi = 1$ and $Wi = 3$ when the steady state is reached. The velocity field exhibits some of the classical features that have been observed in the literature. Each corner at the bottom of the cavity contains a vortex at very low velocity that effectively behaves as a stagnant zone. As we increase the $Wi$, we transition from a left-right symmetry ($Wi = 0$) to an asymmetric profile with a shift of the primary central vortex and of the region of maximum velocity close to the lid towards the left-hand side. The most interesting features are observed for the trace of the conformations tensor. We see a dramatic increase of the maximum value of the trace from $Wi = 1$ to $Wi = 3$, going from approximately $5 \times 10^3$ to more than $10^5$. This maximum is reached at the lid with zones of large stresses both in the middle of the lid and in the top right corner. This suggests that mesh convergence in these regions may be particularly difficult to obtain. We also observe that regions of high stresses have extended further inside the cavity at $Wi = 3$ compared to $Wi = 1$, which may impact convergence even away from the lid.

We then study mesh convergence and compare quantitatively our scheme to results from the literature, starting with the Weissenberg number equal to 1 and comparing to results in [23]. The computation is performed up to a time $t = 10$, with a steady state reached for $t \approx 8$. The first component of the velocity, $u_x$, along the line $x = 0.5$, is presented in Figure 7a. The steady state is almost independent from the mesh and is in close agreement with results in [23]. The second component of the velocity, $u_y$, along the line $y = 0.75$ is plotted in Figure 7b. Convergence is slower and seems to be roughly obtained for the $256 \times 256$ grid, which is also in close agreement with [23].

The most difficult aspect of these computations is to accurately estimate the conformation tensor near the lid. We investigate this issue here using the Mn meshes that are specifically designed for this purpose, with finer meshes in the region close to the lid. Figure 7c shows tr(c) along the line $x = 0.5$ for M1, M2, M3 and M4. On this plot, convergence seems to be achieved for any of the grids. However, we can see in Figure 7d that the profile of tr(c) becomes very steep at the top. Figures 7e and 7f show tr(c) along the lines $y = 0.975$ and $y = 1$ respectively, for the different meshes. At $y = 0.975$, convergence seems almost achieved. The picture is completely different at $y = 1$ where the maximum value of tr(c), obtained close to $x = 0.5$, increases dramatically when refining the mesh. Convergence is very slow as, even for the finest mesh, we see a multiplication of the value of the maximum by a factor of about 1.5 each time the mesh step is divided by $(n+1)/n$ for $n = 1, \dots, 3$.

We now focus on the case $Wi = 3$ and compare with results in [18,23]. When the mesh is coarse, an unsteady behavior is observed that disappears when refining the mesh. The temporal evolution of the kinetic energy,

$$E_{\text{kin}}(t) = \frac{1}{2} \int_{\Omega} |\mathbf{u}(\mathbf{x}, t)|^2 d\mathbf{x},$$

is presented in Figure 8. For the uniform $64 \times 64$ grid, we observe periodic oscillations of the kinetic energy that are related to the appearance of vortices in the upper-right corner of the cavity. However, this phenomenon seems to be purely numerical, since, with finer meshes, oscillations disappear and a stationary solution is obtained after $t \approx 40$. Such numerical instabilities had also been observed in [23], although they did not obtain a stable solution by refining the mesh. In [18], solutions for this case may also be slightly unsteady. We compare our stationary solution to the solutions obtained in these two papers at a given time $t = 40$. Figures 9a and 9b depict profiles of the velocity field for the Mn meshes, showing that convergence is slower than for $Wi = 1$. We find that our results for the first velocity component (Fig. 9a) are in very good agreement with those of
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Evolution of velocity and tr(c) fields with the Weissenberg number. Lid-driven cavity, Oldroyd-B model. From top to bottom: $Wi = 1$, $Wi = 3$ and $Wi = 5$ (with RT discretization). Left: velocity magnitude and streamlines. Right: trace of the conformation tensor (in log-scale). Results obtained with mesh M4.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>FIGURE 6. Evolution of velocity and tr(c) fields with the Weissenberg number. Lid-driven cavity, Oldroyd-B model. From top to bottom: $Wi = 1$, $Wi = 3$ and $Wi = 5$ (with RT discretization). Left: velocity magnitude and streamlines. Right: trace of the conformation tensor (in log-scale). Results obtained with mesh M4.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

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

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

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

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

<figure><img src="image_13.png" /><figcaption>FIGURE 7. Profiles of $u_x$, $u_y$ and $tr(c)$ along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 1$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_8.png" /></figure>

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

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

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

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

<figure><img src="image_13.png" /><figcaption>FIGURE 7. Profiles of $u_x$, $u_y$ and $tr(c)$ along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 1$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIGURE 8. Time evolution of the kinetic energy $E_{\text{kin}}$ for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 3$. (a) Uniform grids. (b) $Mn$ meshes (insert: zoom between $t = 40$ and $t = 80$).</figcaption></figure>

the literature [18, 23]. For the second component, we observe differences with both [18, 23] in the region where $x \gtrsim 0.5$. However, both results in [18, 23] are also significantly different from one another and our solution lies roughly in between the results of these two papers. The zone where this discrepancy is observed is also the area which is likely to be impacted by unstationarities generated near the upper right corner (this was observed in our computations with coarse meshes). For the conformation tensor, difficulties in converging are no longer confined to the very top of the cavity but rather propagate further inwards (see Figs. 9c–9e). As for $Wi = 1$, we still obtain an acceptable convergence everywhere except in the region closest to the lid (Fig. 9f).

Since momentum and conformation transport are strongly coupled, the absence of convergence for the conformation tensor locally in the lid region, which we believe is a common feature of all numerical results so far, raises concern for the convergence of the velocity field. The first clue that we may have convergence for the velocity field is of course the observed mesh convergence for the velocity components, together with the
<figure><img src="image_14.png" /><figcaption>FIGURE 8. Time evolution of the kinetic energy $E_{\text{kin}}$ for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 3$. (a) Uniform grids. (b) $Mn$ meshes (insert: zoom between $t = 40$ and $t = 80$).</figcaption></figure>

the literature [18, 23]. For the second component, we observe differences with both [18, 23] in the region where $x \gtrsim 0.5$. However, both results in [18, 23] are also significantly different from one another and our solution lies roughly in between the results of these two papers. The zone where this discrepancy is observed is also the area which is likely to be impacted by unstationarities generated near the upper right corner (this was observed in our computations with coarse meshes). For the conformation tensor, difficulties in converging are no longer confined to the very top of the cavity but rather propagate further inwards (see Figs. 9c–9e). As for $Wi = 1$, we still obtain an acceptable convergence everywhere except in the region closest to the lid (Fig. 9f).

Since momentum and conformation transport are strongly coupled, the absence of convergence for the conformation tensor locally in the lid region, which we believe is a common feature of all numerical results so far, raises concern for the convergence of the velocity field. The first clue that we may have convergence for the velocity field is of course the observed mesh convergence for the velocity components, together with the
<figure><img src="image_15.png" /><figcaption>FIGURE 9. Profiles of $u_x$, $u_y$ and $tr(c)$ along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 3$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>FIGURE 9. Profiles of $u_x$, $u_y$ and $tr(c)$ along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 3$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
perfect consistency between the MAC and RT results (not shown here). We will also see in the next section that studying a limit of the FENE-CR model provides another clue that the velocity profiles are actually correct.

### 5.1.2. Oldroyd-B model at $Wi = 5$

Finally, we consider the case $Wi = 5$. This case has been studied in [20, 23], where the authors obtained an unsteady regime. In [23], where they use a uniform grid MAC discretization, they obtained oscillations of the kinetic energy with a frequency that decreases with the mesh size. We obtain a similar result in Figure 10 (top) with our MAC scheme on coarse meshes and a stabilization for the M1 mesh (Fig. 10, bottom). These oscillations arise from the fact that the zones where $\text{tr}(c)$ is high form a very thin filament for $Wi = 5$ that is very difficult to capture numerically. When the mesh is not sufficiently refined, the filament can rupture, generating energy transfer and kinetic energy oscillations. Refining the mesh delays this rupture and decreases the frequency of the kinetic energy oscillations. In contrast with what happens at lower Weissenberg number, these oscillations generate slight perturbations which remain apparent on the profile of the $y$-component of the velocity along the line $y = 0.75$, in the right part of the curve (let us say for $0.7 \le x \le 0.9$); except for the conformation tensor very close to the lid, all the other profiles do not vary with time and match the results of the RT discretization (Fig. 12, see comments below). By switching to the RT discretization, oscillations disappear and we visually obtain a plateau of kinetic energy at $t \approx 50$, Figure 11. However, it seems that oscillations of very small amplitude either subsist in the numerical solution or may redevelop if the time step is relaxed. Whatever the reason for this, this prevents increasing the time step at the end of the computation. For the M4 mesh, the minimum time step is set at $10^{-6}$ and the time step remains close to $10^{-5}$ at the end of the computing time. However, in this case also, we believe that the mechanism leading to temporal oscillations is once again numerical.

The results obtained at $Wi = 5$ at steady state are shown on Figure 12 for the RT discretization. As in the previous cases, convergence is obtained on the velocity and the conformation tensor profiles, except close to the lid, for the conformation tensor.

### 5.1.3. FENE-CR model

We now study the solutions of the FENE-CR model. As we will see, convergence is more easily reached than with the Oldroyd-B model, so the MAC and RT discretizations yield equivalent results. The solutions plotted here are obtained with the MAC scheme.

The 2D plots of the trace of the conformation tensor (log-scale) for $Wi = 1$ and $Wi = 3$ are compared in Figure 13 for different values of $b$. Zones of large polymeric stresses are similar to that of the Oldroyd-B model, with large elongations in the lid region. However, when normalized by $b$, the regions of large stresses are much more diffuse and their size tends to reduce as $b$ increases. To understand this, recall that the primary difference between the Oldroyd-B and FENE-CR models is that the elongation of the polymer is bounded by $b$ in the case of the FENE-CR model, while it can become arbitrarily large in the Oldroyd-B model. In all results below, this means that the maximum value of $\text{tr}(c)$ is lower than $b$. Small values of $b$ can considerably limit the polymer elongation, the polymeric stress and the impact of polymers on flow. This is also clear when considering the limits of the FENE-CR model, with, at least formally, $b \to d$ leading to Stokes flow and $b \to \infty$ to the Oldroyd-B model.

Figures 14 and 15 show the mesh convergence for $b = 1000$ at $Wi = 1$ and $Wi = 3$, respectively. For the components of the velocity, convergence is similar to that of the Oldroyd-B model. For the conformation tensor, the behavior is radically different, especially close to the lid. Convergence is much faster in the case of the FENE-CR model, as expected from the effects of the limitation of the polymer elongation. Figures 16 and 17 show results of the FENE-CR model for different values of $b$ at $Wi = 1$ and $Wi = 3$, respectively. These are further compared with the Stokes and Oldroyd-B models. As expected, results for the FENE-CR model get closer to the Newtonian solution when $b$ decreases and to the Oldroyd-B solution in the limit $b \to +\infty$. Figure 18 shows
perfect consistency between the MAC and RT results (not shown here). We will also see in the next section that studying a limit of the FENE-CR model provides another clue that the velocity profiles are actually correct.

### 5.1.2. Oldroyd-B model at $Wi = 5$

Finally, we consider the case $Wi = 5$. This case has been studied in [20, 23], where the authors obtained an unsteady regime. In [23], where they use a uniform grid MAC discretization, they obtained oscillations of the kinetic energy with a frequency that decreases with the mesh size. We obtain a similar result in Figure 10 (top) with our MAC scheme on coarse meshes and a stabilization for the M1 mesh (Fig. 10, bottom). These oscillations arise from the fact that the zones where $\text{tr}(c)$ is high form a very thin filament for $Wi = 5$ that is very difficult to capture numerically. When the mesh is not sufficiently refined, the filament can rupture, generating energy transfer and kinetic energy oscillations. Refining the mesh delays this rupture and decreases the frequency of the kinetic energy oscillations. In contrast with what happens at lower Weissenberg number, these oscillations generate slight perturbations which remain apparent on the profile of the y-component of the velocity along the line $y = 0.75$, in the right part of the curve (let us say for $0.7 \le x \le 0.9$); except for the conformation tensor very close to the lid, all the other profiles do not vary with time and match the results of the RT discretization (Fig. 12, see comments below). By switching to the RT discretization, oscillations disappear and we visually obtain a plateau of kinetic energy at $t \approx 50$, Figure 11. However, it seems that oscillations of very small amplitude either subsist in the numerical solution or may redevelop if the time step is relaxed. Whatever the reason for this, this prevents increasing the time step at the end of the computation. For the M4 mesh, the minimum time step is set at $10^{-6}$ and the time step remains close to $10^{-5}$ at the end of the computing time. However, in this case also, we believe that the mechanism leading to temporal oscillations is once again numerical.

The results obtained at $Wi = 5$ at steady state are shown on Figure 12 for the RT discretization. As in the previous cases, convergence is obtained on the velocity and the conformation tensor profiles, except close to the lid, for the conformation tensor.

### 5.1.3. FENE-CR model

We now study the solutions of the FENE-CR model. As we will see, convergence is more easily reached than with the Oldroyd-B model, so the MAC and RT discretizations yield equivalent results. The solutions plotted here are obtained with the MAC scheme.

The 2D plots of the trace of the conformation tensor (log-scale) for $Wi = 1$ and $Wi = 3$ are compared in Figure 13 for different values of $b$. Zones of large polymeric stresses are similar to that of the Oldroyd-B model, with large elongations in the lid region. However, when normalized by $b$, the regions of large stresses are much more diffuse and their size tends to reduce as $b$ increases. To understand this, recall that the primary difference between the Oldroyd-B and FENE-CR models is that the elongation of the polymer is bounded by $b$ in the case of the FENE-CR model, while it can become arbitrarily large in the Oldroyd-B model. In all results below, this means that the maximum value of $\text{tr}(c)$ is lower than $b$. Small values of $b$ can considerably limit the polymer elongation, the polymeric stress and the impact of polymers on flow. This is also clear when considering the limits of the FENE-CR model, with, at least formally, $b \to d$ leading to Stokes flow and $b \to \infty$ to the Oldroyd-B model.

Figures 14 and 15 show the mesh convergence for $b = 1000$ at $Wi = 1$ and $Wi = 3$, respectively. For the components of the velocity, convergence is similar to that of the Oldroyd-B model. For the conformation tensor, the behavior is radically different, especially close to the lid. Convergence is much faster in the case of the FENE-CR model, as expected from the effects of the limitation of the polymer elongation. Figures 16 and 17 show results of the FENE-CR model for different values of $b$ at $Wi = 1$ and $Wi = 3$, respectively. These are further compared with the Stokes and Oldroyd-B models. As expected, results for the FENE-CR model get closer to the Newtonian solution when $b$ decreases and to the Oldroyd-B solution in the limit $b \to +\infty$. Figure 18 shows
<figure><img src="image_16.png" /><figcaption>FIGURE 10. Time evolution of the kinetic energy $E_{kin}$ for different meshes with MAC discretization. Lid-driven cavity, Oldroyd-B model, $Wi = 5$. (a) Uniform grids. (b) Mn meshes (insert: zoom between $t = 40$ and $t = 100$).</figcaption></figure>

the evolution, as a function of $b$, of the difference (in L2-norm) between the velocity fields:

$$ \| u_b - u_O \|_{L^2} = \left( \int_{\Omega} \| u_b(x) - u_O(x) \|^2 dx \right)^{1/2}, $$

with $u_O$ the velocity field obtained with the Oldroyd-B model and $u_b$ that obtained with the FENE-CR model (for this norm computation, each velocity component is supposed to be piecewise constant over its associated dual mesh). By increasing $b$, we converge towards the Oldroyd-B model, with a faster convergence at $Wi = 1$. Since results for the FENE-CR model are converged for both the velocity and conformation fields, this convergence supports the claim that the velocity profiles for the Oldroyd-B fluid are correct.

## 5.2. Flow past a confined cylinder

We now consider the case of a two-dimensional flow around a confined cylinder. Schematics of the geometry are presented in Figure 19. It consists of a cylinder of radius $R = 1$ confined in a channel of width $W = 4$ and
<figure><img src="image_16.png" /><figcaption>FIGURE 10. Time evolution of the kinetic energy $E_{kin}$ for different meshes with MAC discretization. Lid-driven cavity, Oldroyd-B model, $Wi = 5$. (a) Uniform grids. (b) Mn meshes (insert: zoom between $t = 40$ and $t = 100$).</figcaption></figure>

the evolution, as a function of $b$, of the difference (in L2-norm) between the velocity fields:

$$ \| u_b - u_O \|_{L^2} = \left( \int_{\Omega} \| u_b(x) - u_O(x) \|^2 dx \right)^{1/2}, $$

with $u_O$ the velocity field obtained with the Oldroyd-B model and $u_b$ that obtained with the FENE-CR model (for this norm computation, each velocity component is supposed to be piecewise constant over its associated dual mesh). By increasing $b$, we converge towards the Oldroyd-B model, with a faster convergence at $Wi = 1$. Since results for the FENE-CR model are converged for both the velocity and conformation fields, this convergence supports the claim that the velocity profiles for the Oldroyd-B fluid are correct.

## 5.2. Flow past a confined cylinder

We now consider the case of a two-dimensional flow around a confined cylinder. Schematics of the geometry are presented in Figure 19. It consists of a cylinder of radius $R = 1$ confined in a channel of width $W = 4$ and
<figure><img src="image_17.png" /><figcaption>FIGURE 11. Time evolution of the kinetic energy $E_{\text{kin}}$ for the Mn meshes with RT discretization. Lid-driven cavity, Oldroyd-B model, $Wi = 5$.</figcaption></figure>

length $L = 30$. No-slip and no-penetration boundary conditions, $u_x = u_y = 0$, are imposed on the top/bottom boundary and on the cylinder. Periodic boundary conditions are used for the left/right boundaries. This yields a top/bottom symmetry, which allows us to simulate the flow only in the upper half of the channel. The flow results from a uniform body force $F = F e_x$ (where $e_x$ stands for the first vector of the canonical basis of $\mathbb{R}^2$), with $F$ adjusted so that the flow rate is constant and set to $Q = 2$. We further fix $\beta = 0.59$ and perform computations for $Wi = 0.6$ and $Wi = 0.7$.

The presence of a stagnation point in the wake of the cylinder generates a thin zone of very large polymer elongation, a “filament”, located on the symmetry axis and starting from the intersection of the latter with the boundary of the cylinder. This filament is often referred to as a birefringent strand in the literature, in reference to an experimental method to visualize it. It is notoriously difficult to capture in flow simulations because of both the large values of the conformation tensor and the strong gradients. Mesh convergence for the conformation tensor $c$ and, equivalently, for the polymer stress tensor, requires a very refined discretization in this area. To solve this issue, we use a set of non-structured meshes (Fig. 19) with the RT space approximation. These meshes are built with the GMSH mesher (version 4.4) as follows:

- outside of the blue zone in Figure 19 (delimited by the vertices $(12.5, 0)$, $(12.5, 0.5)$, $(14, 2)$, $(16, 2)$, $(17.5, 0.5)$ and $(17.5, 0)$, if the bottom left corner of the channel is associated with the origin), the length of the cell edges is roughly kept constant at value $\delta_{\text{max}}$ and cells have two horizontal edges. Mesh elements are almost square, except near the symmetry axis where they take the form of (rather flat) rectangles, since the refinement described in the following item propagates up to the left and right boundaries.

- inside the blue zone, the mesh is refined. Cells always have two edges aligned with the radial direction of the cylinder. Along this direction, the size of the cells progressively increases from $\delta_{\text{min}}$ at the boundary to a value close to $\delta_{\text{max}}$ at the blue line that delimits the refined zone. Along the cylinder boundary, the size of the cell edges progressively increases and then decreases, the smallest cells being located at the intersections between the cylinder boundary and the symmetry axis, and being almost square.

Two families of meshes are built in this way. For the first one, meshes are denoted by MC1, ..., MC5, $\delta_{\text{max}} = 1/100$ and the total number of cells is close to 400 000. For the second one, meshes are denoted by MR1, ..., MR7, $\delta_{\text{max}} = 1/300$ and the total number of cells is close to 4 millions. For each of these meshes, the smallest cell width $\delta_{\text{min}}$ is given in Table 1a (MC family of meshes) and Table 1b (MR family of meshes). The time step is chosen between $10^{-3}$ (for the coarsest meshes) and $3 \times 10^{-5}$ (for the most refined ones) and, as in the case of the lid-driven cavity, the sub-time-step for the solution of the ODE on the cell $K$ is set to $\delta t/n_e$.
<figure><img src="image_17.png" /><figcaption>FIGURE 11. Time evolution of the kinetic energy $E_{\text{kin}}$ for the Mn meshes with RT discretization. Lid-driven cavity, Oldroyd-B model, $W_i = 5$.</figcaption></figure>

length $L = 30$. No-slip and no-penetration boundary conditions, $u_x = u_y = 0$, are imposed on the top/bottom boundary and on the cylinder. Periodic boundary conditions are used for the left/right boundaries. This yields a top/bottom symmetry, which allows us to simulate the flow only in the upper half of the channel. The flow results from a uniform body force $F = F e_x$ (where $e_x$ stands for the first vector of the canonical basis of $\mathbb{R}^2$), with $F$ adjusted so that the flow rate is constant and set to $Q = 2$. We further fix $\beta = 0.59$ and perform computations for $W_i = 0.6$ and $W_i = 0.7$.

The presence of a stagnation point in the wake of the cylinder generates a thin zone of very large polymer elongation, a “filament”, located on the symmetry axis and starting from the intersection of the latter with the boundary of the cylinder. This filament is often referred to as a birefringent strand in the literature, in reference to an experimental method to visualize it. It is notoriously difficult to capture in flow simulations because of both the large values of the conformation tensor and the strong gradients. Mesh convergence for the conformation tensor $c$ and, equivalently, for the polymer stress tensor, requires a very refined discretization in this area. To solve this issue, we use a set of non-structured meshes (Fig. 19) with the RT space approximation. These meshes are built with the GMSH mesher (version 4.4) as follows:

- outside of the blue zone in Figure 19 (delimited by the vertices $(12.5, 0)$, $(12.5, 0.5)$, $(14, 2)$, $(16, 2)$, $(17.5, 0.5)$ and $(17.5, 0)$, if the bottom left corner of the channel is associated with the origin), the length of the cell edges is roughly kept constant at value $\delta_{\text{max}}$ and cells have two horizontal edges. Mesh elements are almost square, except near the symmetry axis where they take the form of (rather flat) rectangles, since the refinement described in the following item propagates up to the left and right boundaries.

- inside the blue zone, the mesh is refined. Cells always have two edges aligned with the radial direction of the cylinder. Along this direction, the size of the cells progressively increases from $\delta_{\text{min}}$ at the boundary to a value close to $\delta_{\text{max}}$ at the blue line that delimits the refined zone. Along the cylinder boundary, the size of the cell edges progressively increases and then decreases, the smallest cells being located at the intersections between the cylinder boundary and the symmetry axis, and being almost square.

Two families of meshes are built in this way. For the first one, meshes are denoted by MC1, ..., MC5, $\delta_{\text{max}} = 1/100$ and the total number of cells is close to 400 000. For the second one, meshes are denoted by MR1, ..., MR7, $\delta_{\text{max}} = 1/300$ and the total number of cells is close to 4 millions. For each of these meshes, the smallest cell width $\delta_{\text{min}}$ is given in Table 1a (MC family of meshes) and Table 1b (MR family of meshes). The time step is chosen between $10^{-3}$ (for the coarsest meshes) and $3 \times 10^{-5}$ (for the most refined ones) and, as in the case of the lid-driven cavity, the sub-time-step for the solution of the ODE on the cell $K$ is set to $\delta t/n_e$.
<figure><img src="image_18.png" /><figcaption>FIGURE 12. Profiles of $u_x$, $u_y$ and $\operatorname{tr}(c)$ along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 5$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $\operatorname{tr}(c)$ along $x = 0.5$. (d) Profile of $\operatorname{tr}(c)$ along $x = 0.5$ near the top. (e) Profile of $\operatorname{tr}(c)$ along $y = 0.975$. (f) Profile of $\operatorname{tr}(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE 12. Profiles of $u_x$, $u_y$ and tr(c) along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, $Wi = 5$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of tr(c) along $x = 0.5$. (d) Profile of tr(c) along $x = 0.5$ near the top. (e) Profile of tr(c) along $y = 0.975$. (f) Profile of tr(c) along $y = 1$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIGURE 13. Trace of the conformation tensor for $Wi = 1$ (left) and $Wi = 3$ at $t = 40$ (right), in log-scale. Lid-driven cavity, FENE-CR model. From top to bottom: $b = 10$, $b = 100$ and $b = 1000$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIGURE 13. Trace of the conformation tensor for $Wi = 1$ (left) and $Wi = 3$ at $t = 40$ (right), in log-scale. Lid-driven cavity, FENE-CR model. From top to bottom: $b = 10$, $b = 100$ and $b = 1000$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIGURE 14. Profiles of $u_x$, $u_y$ and $\operatorname{tr}(c)$ along different lines for different meshes. Lid-driven cavity, FENE-CR model with $b = 1000$, $W_i = 1$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$. (d) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$ near the top. (e) Profile of $\operatorname{tr}(c)/b$ along $y = 0.975$. (f) Profile of $\operatorname{tr}(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIGURE 14. Profiles of $u_x$, $u_y$ and $\operatorname{tr}(c)$ along different lines for different meshes. Lid-driven cavity, FENE-CR model with $b = 1000$, $W_i = 1$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$. (d) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$ near the top. (e) Profile of $\operatorname{tr}(c)/b$ along $y = 0.975$. (f) Profile of $\operatorname{tr}(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIGURE 15. Profiles of $u_x$, $u_y$ and $\operatorname{tr}(c)$ along different lines for different meshes. Lid-driven cavity, FENE-CR model with $b = 1000$, $W_i = 3$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$. (d) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$ near the top. (e) Profile of $\operatorname{tr}(c)/b$ along $y = 0.975$. (f) Profile of $\operatorname{tr}(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>FIGURE 15. Profiles of $u_x$, $u_y$ and $\operatorname{tr}(c)$ along different lines for different meshes. Lid-driven cavity, FENE-CR model with $b = 1000$, $W_i = 3$. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$. (d) Profile of $\operatorname{tr}(c)/b$ along $x = 0.5$ near the top. (e) Profile of $\operatorname{tr}(c)/b$ along $y = 0.975$. (f) Profile of $\operatorname{tr}(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIGURE 16. Comparison with Oldroyd-B model for different values of the elongation parameter *b*. Lid-driven cavity, FENE-CR model, $Wi = 1$. Results obtained with mesh M4. (a) First velocity component along *x* = 0.5. (b) Second velocity component along *y* = 0.75. (c) Profile of tr(**c**) along *x* = 0.5. (d) Profile of tr(**c**) along *x* = 0.5 near the top. (e) Profile of tr(**c**) along *y* = 0.975. (f) Profile of tr(**c**) along *y* = 1.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>FIGURE 16. Comparison with Oldroyd-B model for different values of the elongation parameter b. Lid-driven cavity, FENE-CR model, $Wi = 1$. Results obtained with mesh M4. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIGURE 17. Comparison with Oldroyd-B model for different values of the elongation parameter $b$. Lid-driven cavity, FENE-CR model, $Wi = 3$. Results obtained with the mesh M4. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIGURE 17. Comparison with Oldroyd-B model for different values of the elongation parameter $b$. Lid-driven cavity, FENE-CR model, $Wi = 3$. Results obtained with the mesh M4. (a) First velocity component along $x = 0.5$. (b) Second velocity component along $y = 0.75$. (c) Profile of $tr(c)$ along $x = 0.5$. (d) Profile of $tr(c)$ along $x = 0.5$ near the top. (e) Profile of $tr(c)$ along $y = 0.975$. (f) Profile of $tr(c)$ along $y = 1$.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIGURE 18. Difference from the Oldroyd-B velocity field for different values of the elongation parameter *b*. Lid-driven cavity, FENE-CR model, $Wi = 1$ and $Wi = 3$. Results obtained with the mesh M4.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 19. Flow past a confined cylinder – Geometry outline and coarse version of the mesh.</figcaption></figure>

<table><caption>TABLE 1. Minimum width $\delta_{min}$ of the cells located in the wake of the cylinder.</caption><thead><tr><th colspan="8" style="text-align:center;">(A) MC-type meshes</th></tr><tr><th colspan="8" style="text-align:center;">Minimum cell size of meshes of type MC</th></tr><tr><th>Mesh</th><th>MC1</th><th>MC2</th><th>MC3</th><th>MC4</th><th>MC5</th><th colspan="2"></th></tr></thead><tbody><tr><td>δ<sub>min</sub></td><td>8.0 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-4</sup></td><td>1.0 × 10<sup>-4</sup></td><td>5.0 × 10<sup>-5</sup></td><td colspan="2"></td></tr></tbody><thead><tr><th colspan="8" style="text-align:center;">(B) MR-type meshes</th></tr><tr><th colspan="8" style="text-align:center;">Minimum cell size of meshes of type MR</th></tr><tr><th>Mesh</th><th>MR1</th><th>MR2</th><th>MR3</th><th>MR4</th><th>MR5</th><th>MR6</th><th>MR7</th></tr></thead><tbody><tr><td>δ<sub>min</sub></td><td>6.0 × 10<sup>-4</sup></td><td>3.0 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-4</sup></td><td>8.0 × 10<sup>-5</sup></td><td>4.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-5</sup></td><td>1.0 × 10<sup>-5</sup></td></tr></tbody></table>

with $n_e$ the smallest integer number such that $\delta t/n_e \le 1/(2m ||\nabla_K u^{n+1}||_\infty)$ for the Oldroyd-B model and $\delta t/n_e \le (b-2)/(4\lambda b m ||\nabla_K u^{n+1}||_\infty^2)$ for the FENE-CR model, with $m = 100$.
<figure><img src="image_24.png" /><figcaption>FIGURE 18. Difference from the Oldroyd-B velocity field for different values of the elongation parameter b. Lid-driven cavity, FENE-CR model, $Wi = 1$ and $Wi = 3$. Results obtained with the mesh M4.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>FIGURE 19. Flow past a confined cylinder – Geometry outline and coarse version of the mesh.</figcaption></figure>

<table><caption>TABLE 1. Minimum width $\delta_{min}$ of the cells located in the wake of the cylinder.</caption><thead><tr><th colspan="8" style="text-align:center;">(A) MC-type meshes</th></tr><tr><th colspan="8" style="text-align:center;">Minimum cell size of meshes of type MC</th></tr><tr><th>Mesh</th><th>MC1</th><th>MC2</th><th>MC3</th><th>MC4</th><th>MC5</th><th></th><th></th></tr></thead><tbody><tr><td>δ<sub>min</sub></td><td>8.0 × 10<sup>-4</sup></td><td>4.0 × 10<sup>-4</sup></td><td>2.0 × 10<sup>-4</sup></td><td>1.0 × 10<sup>-4</sup></td><td>5.0 × 10<sup>-5</sup></td><td></td><td></td></tr></tbody><thead><tr><th colspan="8" style="text-align:center;">(B) MR-type meshes</th></tr><tr><th colspan="8" style="text-align:center;">Minimum cell size of meshes of type MR</th></tr><tr><th>Mesh</th><th>MR1</th><th>MR2</th><th>MR3</th><th>MR4</th><th>MR5</th><th>MR6</th><th>MR7</th></tr></thead><tbody><tr><td>δ<sub>min</sub></td><td>6.0 × 10<sup>-4</sup></td><td>3.0 × 10<sup>-4</sup></td><td>1.5 × 10<sup>-4</sup></td><td>8.0 × 10<sup>-5</sup></td><td>4.0 × 10<sup>-5</sup></td><td>2.0 × 10<sup>-5</sup></td><td>1.0 × 10<sup>-5</sup></td></tr></tbody></table>

with $n_e$ the smallest integer number such that $\delta t/n_e \le 1/(2m ||\nabla_K u^{n+1}||_\infty)$ for the Oldroyd-B model and $\delta t/n_e \le (b-2)/(4\lambda b m ||\nabla_K u^{n+1}||_\infty^2)$ for the FENE-CR model, with $m = 100$.
<figure><img src="image_26.png" /><figcaption>FIGURE 20. Velocity magnitude, pressure and $\operatorname{tr}(\boldsymbol{c})$ (in log-scale) fields. Flow past a confined cylinder, Oldroyd-B model, $\operatorname{Wi} = 0.7$. Top to bottom: velocity magnitude, pressure and trace of the conformation tensor (in log-scale).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 21. Profile of first normal polymer stress $\tau_{xx}$ along the curve coordinate $s$. Flow past a confined cylinder, Oldroyd-B model, $\operatorname{Wi} = 0.6$. Left: profile along the cylinder surface and in the wake of the cylinder. Right: the same profile zoomed in the wake of the cylinder. (a) Profile along the curve coordinate. (b) Profile zoomed in the wake of the cylinder.</figcaption></figure>

### 5.2.1. Oldroyd-B model

We present the results for two Weissenberg numbers, $\operatorname{Wi} = 0.6$ and $\operatorname{Wi} = 0.7$. In both cases, the overall characteristics of the flow are the same so, for two-dimensional plots, we restrict ourselves to the two-dimensional fields (velocity, pressure and trace of the conformation tensor) on the MC4 mesh for the highest Weissenberg number, once the flow has reached the stationary state (Fig. 20). For $\operatorname{Wi} = 0.7$, this occurs at time $t \approx 3$. We observe three regions of large polymer stresses. As discussed previously, a filament is present in the wake of the cylinder. Two large elongation zones also appear between the upper part of the cylinder and the opposite solid surface, at the neighborhood of the boundary, on both sides.