Flows

O Mokhtari
email: omar.mokhtari@toulouse-inp.fr

Y Davit

J.-C Latché

M Quintard
email: michel.quintard@imft.fr

A STAGGERED PROJECTION SCHEME FOR VISCOELASTIC

Keywords: 2010 AMS Subject Classification. 35Q31, 65N12, 65M12, 76M10, 76M12

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 quandrangular 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.

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

ρ(∂ t u + ξu • ∇u) = -∇p + div τ s + div τ p , in (0, T ) × Ω, (1a)
with τ s = η s (∇u + (∇u) t ) and τ p = η p λ f (c)(c -I d ), div u = 0, in (0, T ) × Ω, (1b)

∂ t c + u • ∇c = (∇u) c + c (∇u) t - 1 λ g(c)(c -I d ), in (0, T ) × Ω, ( 1c
)
where Ω is an open bounded connected subset of R d with d ∈ {2, 3}, (0, T ) a finite time interval with T > 0, u the velocity field of the fluid, p the pressure, and c the symmetric positive-definite polymer conformation tensor. The constant coefficients η s , η p , ρ, and λ are respectively the solvent viscosity, the polymer viscosity, the fluid density and the polymer retardation time. The coefficient ξ serves to artificially eliminate u • ∇u, so as to obtain the unsteady Stokes equations for ξ = 0 and the Navier-Stokes equations for ξ = 1. The tensor τ s is the part of the stress accounting for the fluid solvent and τ p is the part of the stress accounting for the presence of polymers. The functions f (c) and g(c) vary depending on the viscoelastic constitutive equations (see [START_REF] Renardy | Mathematical analysis of viscoelastic flows[END_REF] and [START_REF] Bird | Constitutive equations for polymeric liquids[END_REF] for a review) with, for example, the Oldroyd-B model given by f (c) = g(c) = I d and the FENE-CR model by f (c) = g(c) = b b-tr(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 = ρ U L/(η s + η p ), and the Weissenberg number, W i = λU/L, where U and L are a characteristic velocity and lengthscale. 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 [START_REF] Alves | Numerical methods for viscoelastic fluid flows[END_REF] 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 limts of the HWNP was the log-formulation proposed in [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] and used in many works [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF][START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF][START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF][START_REF] Saramito | On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids[END_REF][START_REF] Tomé | Application of the log-conformation tensor to threedimensional time-dependent free surface flows[END_REF]. 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 [START_REF] Balci | Symmetric factorization of the conformation tensor in viscoelastic fluid models[END_REF] or the kernel conformation transformation [START_REF] Afonso | The kernel-conformation constitutive laws[END_REF]. Despite these advances, the HWNP is also often limited 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 ) × (0, T ),

∂ t c + ∂ x c = a(x)c, c (x
, 0) = c 0 (x), x ∈ (0, L 1 ) and c(0, t) = c i , t ∈ (0, T ),

where a(x) is here to mimick the velocity gradients. This problem is similar to that of [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] 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 [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF][START_REF] Leveque | A Linear Hyperbolic System with Stiff Source Terms[END_REF][START_REF] Lorin | Convection systems with stiff source terms[END_REF] that are found, for example, in modelling shallow-water flows with non-flat bottoms [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for Shallow Water flows[END_REF][START_REF] Duran | Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes[END_REF]. A classical method, such as a centered discretization of the source term, yields non-physical solutions [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]. 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 [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF][START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF].

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) [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF][START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF][START_REF] Bouchut | Upwinding of the source term at interfaces for Euler equations with high friction[END_REF][START_REF] Th | Upwinding sources at interfaces in conservation laws[END_REF]. 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 liddriven cavity problem, at least at moderate Wi numbers) may be schematically understood by considering the steady transport problem ∂ 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 [START_REF] Lee | New formulations, positivity preserving discretizations and stability analysis for non-Newtonian flow models[END_REF][START_REF] Lozinski | An energy estimate for the Oldroyd-B model: theory and applications[END_REF] 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 [START_REF] Afonso | The kernel-conformation constitutive laws[END_REF], 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, the Helmholtz free-energy must therefore decrease in time (see [START_REF] Hu | New entropy estimates for the Oldroyd-B model and related models[END_REF] and [START_REF] Barrett | Existence of global weak solutions for some polymeric flow models[END_REF][START_REF] Barrett | Existence of global weak solutions to some regularized kinetic models for dilute polymers[END_REF]) and its time derivative should be negative. An ideal numerical scheme should preserve this inequality at the discrete level and, in so doing, ensure the long-time stability of solutions. Following this idea, [START_REF] Barrett | Existence and approximation of a (regularized) Oldroyd-B model[END_REF][START_REF] Barrett | Finite element approximation of the FENE-P model[END_REF][START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] developed numerical schemes that verify the free-energy estimate and proved the global-in-time existence of solutions to those discretisations, 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 [START_REF]Lid-driven-cavity simulations of Oldroyd-B models using free-energy-dissipative schemes[END_REF] (see also [START_REF] Lukáčová-Medvid'ová | Energy dissipative characteristic schemes for the diffusive Oldroyd-B viscoelastic fluid[END_REF]) that this approach is not sufficient to guarantee accuracy.

Here, we aim for a numerical scheme that (1) satisfies a discrete equivalent of the the free-energy estimate, (2) minimizes CPU costs and is well suited to high-performance computations, (3) shows a good accuracy, in particular to capture steady-state solutions, (4) can be readily used for a variety of viscoelastic constitutive laws.

We present in Sections 2 and 3 an algorithm that satisfies point [START_REF] Afonso | The kernel-conformation constitutive laws[END_REF]. We use a piecewise constant space approximation for the conformation tensor, which offers greats advantages from an algorithmic point of view. This discretization makes it possible to treat the right-hand side of the transport-reaction equation for the conformation tensor (1c) separately and yields a collection of decoupled cell-by-cell problems. 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 [START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], or a nonconforming low-order finite element, namely the Rannacher and Turek (RT) element [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF] for quadrilateral or hexahedric cells (extension to simplicial cells with the Crouzeix-Raviart element [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] is straightforward). To keep the algorithm efficent and robust, we introduce a fractional-step algorithm where the solution of the Navier-Stokes equations (1a)-(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 div τ 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 fromulation enjoys a built-in Lax-Wendroff weak consistency property [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF].

Unfortunately, respecting Points (2) and (3) 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 consistency for the free-energy estimate. Consequently, Point (1) 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 [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], however, this transformation is only partial. In the spirit of [START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF], 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 [START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF], 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, so that it can be easily modified to capture a variety of viscoelastic behaviors (Point (4)). To further gain in efficiency, we use a local time step for the solution of the ODE on each mesh element, thus ensuring stability and preventing a blowup of the CPU cost. 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 only appear 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 3 S platform developed at the French Institut de Radioprotection et Sûreté Nucléaire (IRSN) [START_REF][END_REF] 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. This also allows us to discuss the convergence of the FENE-CR model towards the Oldroyd-B model when the parameter b takes large values and the regularizing effects brought by the FENE-CR model.

Meshes and unknowns

Let M be a decomposition of the domain Ω, supposed to be regular in the usual sense of the finite element literature. The cells may be:

-convex quadrilaterals (d = 2) or hexahedra (d = 3) for a general domain Ω.

-rectangles (d = 2) or rectangular parallelepipeds (d = 3) for a domain Ω with boundaries that are hyperplanes normal to a coordinate axis. E and E(K) are, respectively, the set of all (d -1)-faces σ of the mesh and of the element K ∈ M. The set of faces included in the boundary of Ω is E ext and the set of internal faces is

E int = E \ E ext . A face σ ∈ E int separating cells K and L is denoted by σ = K|L. n K,σ is the outward normal vector to a face σ of K . For 1 ≤ i ≤ d, we denote by E (i) , E (i) int and E (i)
ext the subset of the faces of E, E int and E ext , respectively, that are perpendicular to the i th unit vector of the canonical basis of R d . For K ∈ M and σ ∈ E, |K| is the measure of K and |σ| the (d -1)-measure of the face σ.

The space discretization is staggered, using either the Marker-And Cell (MAC) scheme [START_REF] Harlow | A numerical fluid dynamics calculation method for all flow speeds[END_REF][START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF] or nonconforming low-order finite element approximations, namely the Rannacher and Turek (RT) element [START_REF] Rannacher | Simple nonconforming quadrilateral Stokes element[END_REF] for quadrilateral or hexahedral meshes. The degrees of freedom for the pressure and the conformation (i.e. the discrete pressure and conformation unknowns),

p K , c K , K ∈ M ,
are associated with the cells of the mesh 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. The Dirichlet boundary conditions are taken into account by setting the velocity unknowns associated with an external face to zero, so the set of discrete velocity unknowns reads

D σ D σ ′ σ ′ = K |M K L M |σ | σ = K |L ε = D σ |D σ ′
{u σ,i , σ ∈ E int , 1 ≤ i ≤ d}.
-For the MAC discretization, the degrees of freedom for the i th component of the velocity are located at the centre of the faces σ ∈ E (i) , so the whole set of discrete velocity unknowns reads:

u σ,i , σ ∈ E (i) int , 1 ≤ i ≤ 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 ∈ M is a simplex, a rectangle or a cuboid, we define D K,σ for σ ∈ E(K) as the cone with basis σ and with vertex the mass center of K (see Figure 1). We thus obtain a partition of K in m sub-volumes, where m is the number of faces of the cell, each sub-volume having the same measure |D K,σ | = |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,σ is referred to as the half-diamond cell associated with K and σ. For σ ∈ E int , σ = K|L, we now define the diamond cell D σ associated with σ by D σ = D K,σ ∪ D L,σ . -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 ∈ M and σ ∈ E(K), is now the rectangle or rectangular parallelepiped of basis σ and of measure |D K,σ | = |K|/2 (see Figure 2). We denote by |D σ | the measure (area of volume) of the dual cell D σ and by ε = D σ |D σ ′ the face separating two diamond cells D σ and D σ ′ . The set of the (dual) faces of D σ is denoted by Ē(D σ ).

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 E (i) S associated with the degrees of freedom of the i th component of the :

K x i-1 2 ,j x i-3 2 x i-1 2 x i+ 1 2 x i-1 x i y j-3 2 y j-1 2 y j y j+ 1 2 y j+ 3 2 u x i-1 2 ,j u x i-3 2 ,j u x i+ 1 2 ,j u x i-1 2 ,j-1 u x i-1 2 ,j+1 : ε x i,j : ε x i-1 4 ,j+ 1 2 h x i-1 2 h x i h y j h y j+ 1 2 Figure 2.
Unknown and dual cell for the x-component of the velocity, notations for staggered discretizations.

velocity (S stands for "scheme"):

E (i) S = E (i)
int for the MAC scheme, E int for the RT scheme.

A free-energy preserving pressure correction scheme

Let us first recall the formal derivation of the energy estimate given in [START_REF] Hu | New entropy estimates for the Oldroyd-B model and related models[END_REF] for a class of viscoelastic models including the Oldroyd-B model. The starting point is System (1) where we set f (c) = g(c) = I d which is the definition of the Oldroyd-B model. For the sake of simplicity and without loss of generality, we assume here that ξ = 1. We consider that all the fields are given initially and that boundary conditions are all no-slip conditions. The idea is that there is no source of energy in the system, either through volume forces or boundary conditions, so that, physically, we expect a constant decrease of the free-energy due to dissipation processes. The first step to showing this consists in taking the inner product of the momentum balance Equation (1a) with u,

∂ t E c + u • ∇E c = -∇p + div τ s (u) + div η p λ c • u, (3)
where E c is the kinetic energy defined as

E c = 1 2 ρ |u| 2 .
Let us also consider the trace of Equation (1c) multiplied by η p /(2 λ). In doing so, we suppose that any tensor c solution of Equation (1c) is symmetric definite positive. We have the relationships tr(∇u c) = ∇u : c and tr(c (∇u) t ) = ∇u : c, thanks to the fact that tr(c) = tr(c t ) and tr(c d) = c : d t for any symmetric matrix c of R d×d and any matrix d of R d×d . We obtain

η p 2 λ ∂ t tr(c) + u • ∇tr(c) = η p λ c : ∇u - η p 2 λ 2 tr(c -I d ). ( 4
)
Since c is supposed to be invertible, we may also multiply Equation (1c) by c -1 to get

c -1 ∂ t c + c -1 u • ∇c = c -1 (∇u) c + (∇u) t - 1 λ (I d -c -1 ). (5)
We are now going to consider the trace of this equation. For the first term on the left-hand side, we know (see [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF]Lemma 1.2]) that tr(c -1 ∂ t c) = ∂ t tr(ln(c)). For the second term on the left-hand side, consider that

u • ∇c = d k=1 u k ∂ k c,
with ∂ k c the tensor obtained by deriving each entry of c with respect to the k th space coordinate. By the same identity as for the first term, we therefore have tr(c -1 u • ∇c) = u • ∇tr(ln(c)). For the first term on the right-hand side, we remark that tr(c -1 ∇u c) = tr(c c -1 ∇u) = tr(∇u) = div(u) = 0 since tr(c d) = tr(d c) for any matrix c and d of R d×d . For the second term on the right-hand side, we also have tr((∇u) t ) = 0. Therefore, taking the trace of Equation ( 5) and multiplying by η p 2λ , we get

η p 2 λ ∂ t tr(ln(c)) + u • ∇tr(ln(c)) = - η p 2 λ 2 tr(I d -c -1 ), (6)
which can also be written as

η p 2 λ ∂ t tr(ln(c) + I d ) + u • ∇tr(ln(c) + I d ) = - η p 2 λ 2 tr(I d -c -1 ), (7)
considering that η p 2 λ ∂ t tr(I d ) + u • ∇tr(I d ) = 0. ( 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

∂ t E + u • ∇E = -∇p + div τ s (u) + div η p λ c • u + η p λ c : ∇u - η p 2 λ 2 tr(c + c -1 -2I d ),
where E is the free-energy defined by

E = 1 2 ρ |u| 2 + η p λ tr(c -ln(c) -I d ).
Finally, integrating over Ω yields

d dt Ω E dx + Ω τ s (u) : ∇u dx + η p 2 λ 2 Ω tr(c + c -1 -2I d ) dx = 0, (9)
by integration by parts and using the no-slip boundary condition. 9) thus provides a stability estimate for the solutions of (1).

Time semi-discrete scheme

Our aim is now to derive a time semi-discretization of System (1) satisfying a semi-discrete analogue of Equation [START_REF] Bermudez | Upwind methods for hyperbolic conservation laws with source terms[END_REF]. Let a uniform partition 0 = t 0 < t 1 < . . . < t N = T of (0, T ) be given, with a constant time step δt. We consider the following semi-discretization in time of System [START_REF] Afonso | The kernel-conformation constitutive laws[END_REF], which belongs to the class of pressure correction schemes, Prediction step -Solve for ũn+1 :

ρ δt ũn+1 -u n + ξ ρ u n • ∇ũ n+1 -η s ∆ũ n+1 + ∇p n - η p λ div c n = 0. ( 10a
)
Correction step -Solve for p n+1 and u n+1 :

ρ δt (u n+1 -ũn+1 ) + ∇(p n+1 -p n ) - η p λ div (c n+1 -c n ) = 0, (10b) div u n+1 = 0, (10c) 1 δt c n+1 -c n + u n • ∇c n+1 = ∇u n+1 c n+1 + c n+1 (∇u n+1 ) t - 1 λ (c n+1 -I d ). ( 10d
)
This time semi-discretization satisfies the stability estimate, as stated in the following result.

Proposition 3.1. Let (u n , p n , c n ) 0≤n≤N be the solution of System [START_REF] Bird | Constitutive equations for polymeric liquids[END_REF]. For 0 ≤ n ≤ N , let us suppose that c n is positive definite and let the discrete free-energy of the solution be defined by

E n (u n , p n , c n , δt) = ρ 2 |u n | 2 + η p 2λ tr(c n -ln(c n ) -I d ) + δt 2 2ρ |∇p n - η p λ div c n | 2 . ( 11
)
Then, for 0 ≤ n ≤ N -1, we have the following stability estimate

Ω E n+1 dx + ρ 2 Ω |ũ -u n | 2 dx + δt Ω η s |∇ũ n+1 | 2 dx + η p 2λ 2 Ω tr(c n+1 -(c n+1 ) -1 -2 I d ) dx ≤ Ω E n dx. ( 12
)
Consequently, the integral over the domain of the discrete free-energy is non-increasing.

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

ũn+1 • ρ δt ũn+1 -u n + ξ ρ u n • ∇ũ n+1 ≥ ρ 2 δt |ũ n+1 | 2 -|u n | 2 + ρ 2 u n • ∇|ũ n+1 | 2 .
Hence, using the fact that u n is divergence-free to recast the transport term as a convection term, thanks to the no-slip boundary condition,

Ω ũn+1 • ρ δt ũn+1 -u n + ξ ρ u n • ∇ũ n+1 dx ≥ Ω ρ 2 δt |ũ n+1 | 2 -|u n | 2 dx. ( 13
)
Thanks once again to the boundary conditions, we have

- Ω ũn+1 • η s ∆ũ n+1 dx = Ω η s |∇ũ n+1 | 2 dx. ( 14
)
Therefore, taking the inner product of Equation (10a) with ũn+1 and integrating over the computational domain yields

ρ 2 δt Ω |ũ n+1 | 2 -|u n | 2 dx + Ω η s |∇ũ n+1 | 2 dx + Ω ũn+1 • ∇p n - η p λ div c n dx ≤ 0. ( 15
)
Let us now rewrite Equation (10b) as

ρ δt u n+1 + ∇p n+1 - η p λ div c n+1 = ρ δt ũn+1 + ∇p n - η p λ div c n .
Taking the square of the L 2 -norm of both sides and multiplying by δt/(2 ρ) yields

ρ 2 δt Ω |u n+1 | 2 dx - Ω η p λ u n+1 • div c n+1 dx + δt 2 ρ Ω |∇p n+1 - η p λ div c n+1 | 2 dx = ρ 2 δt Ω |ũ n+1 | 2 dx + Ω ũn+1 • ∇p n - η p λ divc n dx + δt 2 ρ Ω |∇p n - η p λ divc n | 2 dx.
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)

Ω ∇p n+1 • u n+1 = - Ω p n+1 div(u n+1 ) = 0. ( 16
)
Adding with [START_REF] Boyd | Convex Optimization[END_REF], we obtain the following discrete kinetic energy inequality

ρ 2 δt Ω |u n+1 | 2 -|u n | 2 dx + Ω η s |∇ũ n+1 | 2 dx - Ω η p λ u n+1 • div c n+1 dx + T n+1 -T n ≤ 0, (17)
with, for k = n and k = n + 1,

T k = δt 2 ρ Ω |∇p k - η p λ div c k | 2 dx.
Let us now multiply (10d) by (c n+1 ) -1 to obtain

1 δt (c n+1 ) -1 c n+1 -c n +(c n+1 ) -1 u n •∇c n+1 = (c n+1 ) -1 ∇u n+1 c n+1 +(∇u n+1 ) t - 1 λ (I d -(c n+1 ) -1 ). ( 18
)
For any symmetric definite positive matrices c and d, we have tr(c

-1 (c -d)) ≤ tr ln(c) -tr ln(d) (see [14, Lemma 2.1]), which yields 1 δt tr ln(c n+1 ) -tr ln(c n ) + u n • ∇tr ln(c n+1 ) ≥ tr 1 δt (c n+1 ) -1 c n+1 -c n + (c n+1 ) -1 u n • ∇c n+1 . ( 19
)
Taking the trace of Relation ( 18) and integrating over Ω, by arguments already used in the continuous setting, we thus get

Ω tr ln(c n+1 ) -tr ln(c n ) dx ≥ - Ω 1 λ tr(I d -(c n+1 ) -1 ) dx.
Taking the trace of Equation (10d), adding the discrete kinetic energy balance [START_REF] Chen | Hyperbolic conservation laws with stiff relaxation terms and entropy[END_REF] and subtracting this relation as in the continuous setting yields the desired estimate, thanks to the duality relation

Ω tr(∇u n+1 c n+1 ) dx = Ω ∇u n+1 : c n+1 dx = - Ω u n+1 • div c n+1 dx. (20)

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 ũn+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , ρ δt ũn+1 σ,i -u n σ,i + ξ ρ div σ (ũ n+1 i u n ) -div σ,i τ s (ũ n+1 ) + ∇ σ,i p n -div σ,i τ n p = 0. ( 21a
)
Correction step -Solve for p n+1 , u n+1 and c n+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , ρ δt (u n+1 σ,i -ũn+1 σ,i ) + ∇ σ,i (p n+1 -p n ) -div σ,i τ n+1 p + div σ,i τ n p = 0, (21b) ∀K ∈ M, div K (u n+1 ) = 0, ( 21c
) ∀K ∈ M, 1 δt (c n+1 K -c n K ) + div K (u n c n+1 ) = ∇ K u n+1 c n+1 K + (c n+1 K ) t (∇ K u n+1 ) t - 1 λ (c n+1 K -I d ). (21d)
In the prediction step and in the velocity correction step, for k = n and k = n+ 1, the tensor τ k p is discretized on the primal cells and is computed as a function of the conformation tensor as

(τ k p ) K = η p λ (c k K -I d ), for K ∈ 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 ∇u c

+ c t (∇u) t instead of ∇u c + c (∇u) t ), as is done in Equation (21d). Once the equation is symmetrized, if c n is symmetric, for 1 ≤ i, j ≤ d, the discrete equation for c n+1 i,j
is strictly the same as the equation for c n+1 j,i . 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 timestep depending on the L ∞ -norm of ∇u n+1 . Finally, for this equation and this convection scheme, we need a discrete analogue of [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF] 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 [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF]Section 4.2]. The discrete kinetic energy balance (i.e. the discrete equivalent of ( 13)) is established in [START_REF] Latché | A convergent staggered scheme for the variable density incompressible Navier-Stokes equations[END_REF]Lemma 7.2]. For the MAC scheme, a discretization of the convection term is proposed in [24, Section 2, Equation ( 29)] and its L 2 -stability (i.e., as before, the discrete equivalent of ( 13)) is shown in [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF]Lemma 3.6]. The convection operator takes the following general form

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , ρ div σ (ũ n+1 i u n ) = 1 |D σ | ε∈ Ē(Dσ ) F σ,ε (ũ n+1 i ) ε ,
where F σ,ε is an approximation of the mass flux through ε computed from the fluxes through the primal faces (i.e., for K ∈ M and σ ∈ E(K), the quantities F K,σ = ρ |σ| u K,σ with u K,σ given by (25a) and [START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF] for the RT and MAC schemes respectively) in such a way that a divergence-free constraint is also satisfied on the dual cells D σ (i.e. ε∈ Ē(Dσ ) F σ,ε = 0). The approximation of the i th component of the velocity at the dual face, (ũ n+1 i ) ε , may be chosen centered.

Finally, we need to ensure the duality of the velocity divergence and pressure gradient, together with the duality relations ( 14), ( 16) and [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF]. 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 Section 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 [24, Section 2] and a proof of their duality is given by [START_REF] Gallouët | Convergence of the Marker-And-Cell scheme for the incompressible Navier-Stokes equations on non-uniform grids[END_REF]Lemma 2.4]. For the viscous diffusion term, we propose in [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] a construction that ensures the positivity of the dissipation, i.e. a discrete analogue of Relation [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF]. 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 test 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 div σ,i τ 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.

The discrete total stress divergence term

The aim of this section is to define the divergence term div σ,i (T ) of the total stress tensor T = -pI d + τ s + τ p . We want this quantity to satisfy a duality relation of the form

- d i=1 σ∈E (i) S |D σ | div σ,i T u σ,i = K∈M |K| (T : ∇u) K , (22)
which is the discrete analogue of the identity

- Ω div T • u dx = Ω T : ∇u dx, (23)
or, in other words, 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. The discretization of this divergence term 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.

Unstructured discretization

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

-div σ,i T = 1 |D σ | K∈M K T : ∇ϕ (i) σ dx, (24)
where ϕ

(i)
σ stands for the vector-valued finite element shape function associated with the i th component of the velocity and the face σ. By definition of the RT finite elements, this shape function is ϕ σ e (i) , where ϕ σ is the real-valued function of the approximation space whose mean value is 1 over σ 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 2 1 -x 2 2 } on the reference (0, 1) 2 square (d = 2) and {1, x 1 , x 2 , x 3 , x 2 1 -x 2 2 , x 2 2 -x 2 3 } on the reference (0, 1) 3 cube (d = 3).
By a simple reordering of the sums, [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF] reads

- d i=1 σ∈Eint |D σ | div σ,i T u σ,i = d i=1 σ∈Eint u σ,i K∈M K T : ∇ϕ (i) σ dx = K∈M K T : ∇ d i=1 σ∈Eint u σ,i ϕ (i) σ dx = K∈M K T : ∇u dx = K∈M |K| (T : ∇u) K .
Considering separately the pressure gradient, the diffusion and the elastic extra-stress terms in T , we obtain the following expression for the discrete operators featured in the scheme [START_REF] Duran | Asymptotic preserving scheme for the shallow water equations with source terms on unstructured meshes[END_REF], dropping the time index for

: K i,j ∂ M x u x = u x i+ 1 2 ,j -u x i- ,j h x i : K xy i-1 2 ,j-1 2 ∂ M y u x = u x i-1 2 ,j -u x i- ,j-1 h y j-1 2 x i-3 2 x i-1 2 x i+ 1 2 y j-3 2 y j-1 2 y j+ 1 2 u x i-1 2 ,j u x i+ 1 2 ,j u x i-1 2 ,j-1 h x i h y j-1 2 Figure 3
. Discrete partial derivatives of the x-component of the velocity.

: short

K xy i-1 2 ,j-1 2 ∂ M y u x (x) = u x i-1 2 ,j -u x i-1 2 ,ext h y j-1 2 x i-3 2 x i-1 2 x i+ 1 2 y j-1 2 y j+ 1 2 y j+ 3 2 u x i-1 2 ,j u x i-3 2 ,j u x i+ 1 2 ,j u x i-1 2 ,j+1 h x i h y j-1 2
∀K ∈ M, div K u = 1 |K| K div(u) dx = 1 |K| σ∈E(K) |σ| u K,σ , (25a)
with u K,σ = u σ • n K,σ , ∀σ ∈ E int , σ = K|L, ∇ σ p = 1 |D σ | d i=1 M∈M p M M ∇ϕ (i) σ dx = |σ| |D σ | (p L -p K ) n K,σ , (25b)
For

1 ≤ i ≤ d, ∀σ ∈ E (i) S , div σ,i τ s (ũ) = 1 |D σ | Ω τ s (ũ) : ∇ϕ (i) σ dx, (25c) For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , div σ,i τ p = 1 |D σ | K∈M (τ p ) K : K ∇ϕ (i) σ dx, ( 25d
) ∀K ∈ M, ∇ K u = 1 |K| K ∇u dx. (25e)

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 twodimensions, the discrete analogues of ∂ x u x and ∂ 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 [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF] (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 (25) above. This approach was followed in [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] 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 [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] 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 M = {K i,j , (i, j) ∈ I}, with I ⊂ N 2 with notations for geometrical quantities that are given in Figure 3. 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) ∈ I, let the primal cells be

K i,j = (x i-1 2 , x i+ 1 2 ) × (y j-1 2 , y j+ 1 2
). The discrete derivative involved in the divergence (so, for the velocity x-component, only ∂ M x u x ) is defined over the primal cell by, ∀(i, j) ∈ I and ∀x ∈ K i,j ,

∂ M x u x (x) = u x i+ 1 2 ,j -u x i-1 2 ,j h x i . (26)
-For the other derivatives (so, for the velocity x-component, only ∂ M y 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 xy i-1 2 ,j-1 2 = (x i-1 , x i ) × (y j-1 , y j ). We define I xy the set of pairs (i, j) such that x i-1 2 ,j-1 2 is an internal vertex of the mesh (i.e. a vertex which does not lie on the boundary). Then, for (i, j) ∈ I xy and ∀x ∈ K xy i-1 2 ,j-

1 2 : ∂ M y u x (x) = u x i-1 2 ,j -u x i-1 2 ,j-1 h y j-1 2 . ( 27
)
The only necessary extension of this definition to cope with boundaries concerns the definition of ∂ M y u x over a half vertex-centered cell associated with a vertex lying on a horizontal boundary (see Figure 4). In this case, we use the usual "fictitious cell trick" in order to apply Relation ( 27): 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

∇ M u = ∂ M x u x ∂ M y u x ∂ M x u y ∂ M y u y , τ M s (u) = η s ∇ M u + (∇ M u) t . ( 28
)
Finite-volume test functions -Let us denote by

I x ⊂ N 2 (resp. I y ⊂ N 2 ) the set of pairs (i, j) such that x i-1 2 ,j (resp. x i,j-1 2 )
is the mass center of a vertical (resp. horizontal) internal face of the mesh. For (i, j) ∈ I x , we denote by φ x,(i-1 2 ,j) the test function associated with the degree of freedom of the x-component of the velocity located at x i-1 2 ,j . This discrete function is defined as

(φ x,(i-1 2 ,j) ) x k-1 2 ,ℓ = δ i k δ j ℓ , ∀(k, ℓ) ∈ I x and (φ x,(i-1 2 ,j) ) y k,ℓ-1 2 = 0, ∀(k, ℓ) ∈ I y .
Its non-zero partial derivatives are ∂ M x φ x,(i-1 2 ,j) and ∂ M y φ x,(i-1 2 ,j) with

∂ M x φ x,(i-1 2 ,j) = 1 h x i-1 over K i-1,j , -1 h x i over K i,j , 0 elsewhere, ∂ M y φ x,(i-1 2 ,j) = 1 h y j-1 2 over K xy i-1 2 ,j-1 2 , -1 h y j+ 1 2 over K xy i-1 2 ,j+ 1 2 , 0 elsewhere.
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

∀(i, j) ∈ I x , (∇p) x i-1 2 ,j = -1 |K x i-1 2 ,j | Ω p ∂ M x φ x,(i-1 2 ,j) dx. ( 29
)
Using the expression of the discrete derivatives of the test functions, we get

(∇p) x i-1 2 ,j = -1 |K x i-1 2 ,j | (h x i-1 h y j ) p i-1,j 1 h x i-1 + (h x i h y j ) p i,j -1 h x i = p i,j -p i-1,j h x i-1 2 ,
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

for 1 ≤ i ≤ d, ∀σ ∈ E (i) S , σ = K|L, ∇ σ,i p = |σ| |D σ | (p L -p K ) (n K,σ ) i ,
which is a definition similar to the RT discretization, Equation (25b). The velocity divergence operator immediately follows from the definition of the velocity discrete partial derivatives

∀K ∈ M, div K (u) = 1 |K| K div(u) dx,
so, since the partial derivatives of the velocity involved in the divergence are constant over each primal cell,

∀(i, j) ∈ I, div Ki,j (u) = u x i+ 1 2 ,j -u x i-1 2 ,j h x i + u y i,j+ 1 2 -u y i,j-1 2 h y j ,
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

∀K ∈ M, div K (u) = 1 |K| σ∈E(K) |σ| u K,σ , with, for σ ∈ E (i) , u K,σ = u σ n K,σ • e (i) . (30)
Finally, the discrete divergence of the viscous stress tensor is defined by the following weak formulation

∀(i, j) ∈ I x , -(divτ s (ũ)) x i-1 2 ,j = 1 |K x i-1 2 ,j | Ω τ M s (ũ) : ∇ M φ x,(i-1 2 ,j) dx, (31)
By a standard reordering of the sums, we get

d i=1 E∈E (i) S |D σ | ũσ,i div σ,i τ s (ũ) = Ω τ M s (ũ) : ∇ M ũ dx, (32)
An important consequence of this relation, with the definition of τ M s (ũ) = η s ∇ M ũ + (∇ M ũ) t , is that, thanks to an elementary property of the tensors contraction, we may replace the second tensor ∇ M ũ in [START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF] by its symmetrization ∇ M ũ + (∇ M ũ) t /2. We thus obtain that the dissipation (i.e. the integral in [START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF]) is non-negative, which is a discrete counterpart of [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF]. In addition, the duality Relation [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF] requires that the discretization of ∇ K u in the right-hand side of the transport equation for the conformation tensor (21d) to be

∇ K u = 1 |K| K ∇ M u dx.
This weak formulation of the diffusion term ( 32) is shown in [START_REF] Grapsas | An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations[END_REF] 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 (τ M s ) x,x and (τ M s ) y,y , the viscosity must be associated with the primal cells, while, for the shear components (τ M s ) x,y and (τ M s ) 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 τ p to a piecewise function over the primary cells by

∀(i, j) ∈ I, ∀x ∈ K i,j , τ p (x) = η p λ f (c i,j )(c i,j -I d ).
Then we set

∀(i, j) ∈ I x , -(divτ p ) x i-1 2 ,j = 1 |K x i-1 2 ,j | Ω τ p : ∇ M φ x,(i-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 x i-1 2 ,j of a discretization of the flux associated with divτ p , i.e. the integral of the first component of τ p n K,σ (see Figure 2).

-On the vertical face ε = ε x i,j , the flux is given by

ε (τ p n ε ) x = 1 h y j (τ p ) xx i,j = |ε x i,j |(τ p ) xx i,j .
-On the horizontal (half

) face ε = ε x i-1 4 ,j+ 1 2
, the flux is given by

ε (τ p n ε ) x = 1 h y j+ 1 2 |K i,j | 4 (τ p ) xy i,j + |K i,j+1 | 4 (τ p ) xy i,j+1 = |ε x i-1 4 ,j+ 1 2 | α (τ p ) xy i,j + (1 -α)(τ p ) xy i,j+1 , with α = h y j 2h y j+ 1 2 , and 1
-α = h y j+1 2h y j+ 1 2 .
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

ε = ε x i-1 4 ,j+ 1 2
, we have

(τ p ) xy ε = α (τ p ) xy i,j + (1 -α)(τ p ) xy i,j+1
, whereas, with a linear interpolation, we would have

(τ p ) xy ε = (1 -α) (τ p ) xy i,j + α (τ p ) xy i,j+1 ,
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 ∇ 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. This result would extend also to non-linear 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 (Ω) d×d of g(c (m) ) (c (m) -I d )). Lemma 3.2. Let (M (m) ) m∈N be sequence of structured discretizations of Ω, with a space step tending to zero. Let (c (m) ) m∈N be a sequence of discrete conformation tensors, with c (m) associated with M (m) . We suppose that (c (m) ) m∈N weakly converges in

L 1 (Ω) d×d to c. Let ϕ ∈ C ∞ c (Ω) d
, and let us denote by ϕ (m) the interpolate of ϕ in the space of velocity discrete functions defined by ϕ

(m) σ,i = ϕ i (x σ ) for σ ∈ E (i) , with x σ the mass center of σ. Then lim m→+∞ d i=1 Ω div σ,i η p λ (c (m) -I d ) ϕ (m) i dx = - Ω η p λ (c -I d ) : ∇ϕ dx.
Proof. By construction, we have

d i=1 Ω div σ,i η p λ (c (m) -I d ) ϕ (m) i dx = Ω η p λ (c (m) -I d ) : ∇ M ϕ (m) dx,
and the results follows thanks to the convergence of ∇ M ϕ (m) (x) to ∇ϕ(x) for almost every x ∈ Ω.

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 div σ,i τ n+1 p -div σ,i τ n p ) 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)-(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

∇u n+1 K c n+1 K + (c n+1 K ) t (∇u n+1 K ) 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 Section 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 non-linearities. 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 [START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF], 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 ũn+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , ρ δt ũn+1 σ,i -u n σ,i + ξ ρ div σ (ũ n+1 i u n ) -div σ,i τ s (ũ n+1 ) +∇ σ,i p n -div σ,i τ n p = 0. ( 33a
)
Correction step -Solve for p n+1 and u n+1 :

For 1 ≤ i ≤ d, ∀σ ∈ E (i) S , ρ δt (u n+1 σ,i -ũn+1 σ,i ) + ∇ σ,i (p n+1 -p n ) = 0, ( 33b
) ∀K ∈ M, div K (u n+1 ) = 0. ( 33c
)
Up to the (explicit) term div σ,i τ n p 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+ 1 3 :

∀K ∈ M, 1 δt/2 log c n+ 1 3 K -log c n K + div K (u n+1 log c n+ 1 3 ) = 0. (34a) ODE -Set c K (t n ) = c n+ 1 3 K
and solve for c n+ 2 3 = c K (t n + δt):

∀K ∈ M, ∂ t c K -(∇ K u n+1 ) c K -c t K ∇ K u n+1 t + g(c K ) λ (c K -I d ) = 0. ( 34b
)
Advection II -Solve for c n+1 :

∀K ∈ M, 1 δt/2 log c n+1 K -log c n+ 2 3 K + div K (u n+1 log c n+1 ) = 0. ( 34c
)
Transport steps are discretized by a standard backward first-order upwind scheme. From the form of these equations, these steps preserve the symetry of the conformation tensor: if c n is symmetric, the same equation is solved for c n+1 i,j and c n+1 j,i , 1 ≤ i, j ≤ d. In addition, since we search for log(c) and then compute c by taking the exponential of the result, these transport steps necessarily yields positive definite conformation tensors. The local ODE (34b) is solved using a first-order semi-implicit Euler scheme, with a local sub-time step that reads

1 δt K (c n+1 K -c n K ) = ∇u c n+1 K + (c n+1 K ) t (∇u) t - g(c n K ) λ (c n+1 -I d ). ( 35
)
We prove in Section A.2 that this relation yields a symmetric positive definite tensor under the constraint [START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF] for the choice of the sub-time step δt K . In addition, for the FENE-CR model, we prove in Section B that, under additional constraints on the sub-time step (see Lemma B.2), the solution of (35) satisfies tr(c n+1 ) < b.

Numerical results

Here we consider the dimensionless version of the system of equations (1) in the vanishing inertia limit (ξ = 0),

∂ t u = -∇p + β div (∇u + (∇u) t ) + 1 -β W i div f (c) (c -I d ) , (36a)
div u = 0, (36b)

∂ t c + u • ∇c -(∇u) c -c (∇u) t + 1 W i g(c) (c -I d ) = 0, (36c)
where β = η s /(η s + η p ) is the viscosity ratio, W i = λU/L, space is non-dimentionalized with the lengthscale L and time with L/U .

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 [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF][START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF][START_REF] Grillet | Modeling of viscoelastic lid driven cavity flow using finite element simulations[END_REF][START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF][START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF][START_REF] Yapici | Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity[END_REF], each of these cases is challenging to solve numerically even at moderate W i numbers. In the case of the lid-driven cavity, it is only with the recent introduction of the log-conformation formulation that results at W i ≥ 1 have been obtained [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. For the flow around the cylinder, a singularity in the wake of the cylinder appears at W i ≃ 0.7 that is very difficult to capture accurately [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF].

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 +∞. To perform these computations, the proposed scheme was implemented in the open-source CALIF 3 S software developed by the Institut de Radioprotection et de Sûreté Nucléaire IRSN [START_REF][END_REF]. 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 using a Block-Jacobi preconditioner. Parallel computations were then carried out on TotalEnergies' group supercomputer Pangea II.

Lid-driven cavity

In this test, the computational domain is the unit-square Ω = (0, 1) 2 . At the initial time, the fluid is at rest (i.e. u(x, 0) = 0 a.e. in Ω) and the conformation tensor is uniformly equal to the identity (c(x, 0) = I d a.e. in Ω). Dirichlet boundary conditions are imposed to the velocity components at each side of the square, with u t = (8 x 2 (1 -x) 2 1 + tanh(8t -4) , 0) (see [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF]) on the lid and u = 0 on the three other sides and at any time. The viscosity ratio is β = 0.5.

We solve this problem for a sequence of successively refined Cartesian meshes. The three coarsest meshes are uniform 64 × 64, 128 × 128 and 256 × 256 grids. The other ones, denoted by Mn, n = 1, . . . , 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 G (0, 0.001/n), (0.5, 0.004/n) ∪ G (0.5, 0.004/n), (1, 0.001/n) ,

-The set of y-coordinates is G (0, 0.001/n), (1, 0.004/n) , where G((a, δx i ), (b, δx f )) is the set of real numbers beginning at a, ending at b, with a step following a geometrical progression from δx i at a to δ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 δ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 δt = 10 -3 , which is the time step when the solution becomes stationary. The sub-time-step for the solution of the ODE is set to δt/n e with n e the smallest integer number such that δt/n e ≤ 1/(2 m ||∇ K u n+1 || ∞ ) for the Oldroyd-B model and δt/n e ≤ (b -2)/(4 λ b m ||∇ K u n+1 || 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 CPUtime 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. 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 5 the 2D representations of the velocity magnitude and the trace of the conformation tensor for W i = 1 and W i = 3. 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 W i, we transition from a left-right symmetry (W i = 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 W i = 1 to W i = 3, going from approximately 5 × 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 W i = 3 compared to W i = 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 [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. The computation is performed up to a time t = 10, with a steady state reached for t ≈ 8. The first component of the velocity, u x , along the line x = 0.5, is presented in Figure (6a). The steady state is almost independent from the mesh and is in close agreement with results in [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. The second component of the velocity, u y , along the line y = 0.75 is plotted in Figure (6b). Convergence is slower and seems to be roughly obtained for the 256 × 256 grid, which is also in close agreement with [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF].

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 (6c) 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 (6d) that the profile of tr(c) becomes very steep at the top. Figures (6e) and (6f) 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, . . . , 3.

We now focus on the case W i = 3 and compare with results in [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. When the mesh is coarse, an unsteady behavior is observed that disappears when refining the mesh. The temporal evolution of the kinetic energy,

E k = 1 2 Ω |u(x, t)| 2 dx,
is presented in Figure 7. For the uniform grid 64 × 64, 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 ≈ 40. Such numerical instabilities had also been observed in [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], although they did not obtain a stable solution by refining the mesh. In [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF], 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 (8a) and (8b) depict profiles of the velocity field for the Mn meshes, showing that convergence is slower than for W i = 1. We find that our results for the first velocity component (Figure (8a)) are in very good agreement with those of the literature [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. For the second component, we observe differences with both [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF] and [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] in the region where x 0.5. However, both results in [START_REF] Comminal | Robust simulations of viscoelastic flows at high Weissenberg numbers with the streamfunction/log-conformation formulation[END_REF] and in [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] 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 Figures (8c)-(8e)). As for W i = 1, we still obtain an acceptable convergence everywhere except in the region closest to the lid (Figure (8f)).

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

Oldroyd-B model at W i = 5

Finally, we consider the case W i = 5. This case has been studied in [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] and [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF], where the authors obtained an unsteady regime. In [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], 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 9 (top) with our MAC scheme on coarse meshes and a stabilization for the M1 mesh. These oscillations arise from the fact that the zones where tr(c) is high form a very thin filament for W i = 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 left part of the curve (let us say for 0.7 ≤ x ≤ 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 presented on Figure 11. By switching to the RT discretization, oscillations disappear and we visually obtain a plateau of kinetic energy at t ≈ 50, Figure 9 (bottom). 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 W i = 5 at steady state are shown on Figure 11 for the RT discretization. For mesh M4, the time step δt in the final plateau is reduced to avoid a possible numerical destabilization. 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.

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 W i = 1 and W i = 3 are compared in Figure 12 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 tr(c) is lower than b (so the plotted ratio tr(c)/b always remains lower than 1). Small values of b can considerably limit 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 b → d leading to Stokes flow and b → ∞ to the Oldroyd-B model.

Figures 13 and14 show the mesh convergence for b = 1000 at W i = 1 and W i = 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. Figure 15 and 16 show results of the FENE-CR model for different values of b at W i = 1 and W i = 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 → +∞. Figure 17 shows the evolution, as a function of b, of the difference (in L2-norm) between the velocity fields:

u b -u O L 2 = Ω u b (x) -u O (x) 2 dx 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 W i = 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. u tr(c)

W i = 1 W i = 3 W i = 5

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 18. It consists of a cylinder of radius R = 1 confined in a channel of width W = 4 and 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 R 2 ), with F adjusted so that the flow rate is constant and set to Q = 2. We further fix β = 0.59 and obtain the target Weissenberg number, W i = λQ/W R, by changing the relaxation time λ.

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

τ p = η p λ f (c) (c -I d )
requires a very refined discretization in this area. To solve this issue, we use a set of non-structured meshes (Figure 18) 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 18 (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 δ 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 δy min at the boundary to a value close to δ 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, δ max = 1/100 and the total number of cells is close to 400, 000. For the second one, meshes are denoted by MR1, ..., MR7, δ max = 1/300 and the total number of cells is close to 4 millions. For each of these meshes, the smallest cell width, δy 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 × 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 is set to δt/n e with n e the smallest integer number such that δt/n e ≤ 1/(2 m ||∇ K u n+1 || ∞ ) for the Oldroyd-B model and

δt/n e ≤ (b -2)/(4 λ b m ||∇ K u n+1 || 2
∞ ) for the FENE-CR model, with m = 100.

Oldroyd-B model

We present the results for two Weissenberg numbers, W i = 0.6 and W i = 0.7. In both cases, the overall characteristics of the flow are the same. So we only plot in Figure 19 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 stationary state. For W i = 0.7, this occurs at time t ≈ 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. Table 1. Minimum width δy min of the cells located in the wake of the cylinder.

To go further, we compare our results with those proposed in [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF]. For W i = 0.6, the stationary regime is obtained at time t ≈ 2. Figure 20 shows the xx-component of the polymer stress,

τ xx = η p λ (c xx -1),
as a function of the curvilinear coordinate s along the line composed of the cylinder boundary and the part of the symmetry axis located in the wake of the cylinder. The two stagnations points at the cylinder boundary correspond to s = 0 for the front one and s = π for the rear one. As expected from the two-dimensional plot, we observe two peaks for τ xx . The first one, more pronounced, is approximately centered at the throat between the cylinder and the solid surface, s ≈ π/2, and the second one is located in the wake close to the cylinder. These results are in very good agreement with both [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF] and [START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF].

We now turn to the case W i = 0.7, for which a lack of convergence on τ xx in the wake of the cylinder is generally reported in the literature (see [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF]). This is observed in Figure 21 (left) with the MC-type meshes. We find a profile close to the one proposed in [START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF] with the MC1 mesh and close to the one proposed in [START_REF] Damanik | A monolithic FEM approach for the log-conformation reformulation (LCR) of viscoelastic flow problems[END_REF] with the MC3 one. Beyond MC3, we need finer meshes that make it possible to capture thinner filaments with larger values of the conformation tensor. To explore this, we consider the second family of refined meshes, MR1, ..., MR7. The stress profiles are shown in Figure 21 (right), making it obvious that mesh convergence is very difficult to obtain and that the peak value increases as the mesh is refined. However, convergence seems close with the MR7 mesh. Note that, to simplify notations, we do not associate a time-index to the velocity, but the fact that the velocity depends on time is indifferent in the present proof. Let X ∈ R d be a (fixed) vector. Then, for K ∈ M and 1 ≤ i, j ≤ d, we have

FENE-CR model

|K| δt X i X j (c i,j ) n+1 K -X i X j (c i,j ) n K + σ=K|L |σ| X i X j (c i,j ) n+1 σ u K,σ = 0.
Let us define S n K = (c n K X, X). Then, summing the previous relation over i and j yields

|K| δt (S n+1 K -S n K ) + σ=K|L |σ| S n+1 σ u K,σ = 0, with, thanks to Definition (39), S n+1 σ = S n+1 K if u K,σ ≥ 0 and S n+1 σ = S n+1
L otherwise. We may see S as an unknown obeying a first-order backward-Euler upwind scheme, which is known to preserve the positivity. For any reasonable initialization, we have S 0 K > 0 for K ∈ M, so S n K > 0 for n > 0, which means that the discrete tensor c is positive in each cell and at any time.

Remark A.2. In the case with the change of variable to the log-conformation tensor, taking the exponential after the advection step allows to automatically ensure the positivity of c.

A.2. The pure reaction problem

With the implicit variant, the ODE step for the conformation tensor can be written in each cell, dropping for short the cell index,

1 δt (c n+1 -c n ) = ∇u c n+1 + (c n+1 ) t (∇u) t - 1 λ (c n+1 -I d ). (40)
Note that, in the second term on the right-hand side, we have replaced c n+1 by (c n+1 ) t . This is justified by the fact that, provided that c n is symmetrical, any solution of ( 40) is also symmetrical and thus satisfies the initial equation. Equation ( 40) is a linear system for the components of c n+1 which is regular for δt small enough. For instance, the condition

1 δt + 1 λ > 2 max 1≤i≤d d j=1 |(∇u) i,j )| (41)
ensures that the diagonal of the matrix of this system is strictly dominant. The following lemma states that the scheme (40) preserves the definite positivity of the unknown.

Lemma A.3. Let us suppose that the time step is such that Condition (41) is satisfied, and that c n is a symmetric positive definite tensor. Then so is c n+1 .

Proof. Let α → c(α), [0, 1] → R d×d be the function defined by taking c(α) as the solution to Equation [START_REF] Lukáčová-Medvid'ová | Energy dissipative characteristic schemes for the diffusive Oldroyd-B viscoelastic fluid[END_REF] with α δt instead of δt:

1 δt c(α) -c n = α ∇u c(α) + c(α) t (∇u) t - 1 λ c(α) -I d . (42)
Since, thanks to Condition (41), c(α) is the solution to a regular linear system, the function c(α) is regular. In addition, c(0) = c n is a symmetric positive definite tensor. Let us suppose that c(α) is not definite positive for some α ∈ [0, 1], which means that one at least of its eigenvalues is non-positive. Then, by continuity of α → min(eigenvalues of c(α)), there exists α 0 ∈ (0, 1] such that c(α 0 ) has a zero eigenvalue. Let X ∈ R d be where b is a real number strictly greater than the space dimension d. From Section A.2, we easily deduce that, provided that c n is a symmetric definite positive tensor such that tr(c n ) < b, then c n+1 is also symmetric definite positive, under the stability constraint for the time step:

1 δt > 2 max 1≤i≤d d j=1 |(∇u) i,j )|. (46)
We search for an upper bound for δt that ensures that tr(c n+1 ) < b. To this purpose, we will use the following lemma.

Lemma B.1. Let A be a symmetric positive d × d matrix of entries (A i,j ) 1≤i,j≤d . Then:

||A|| L 1 = d i,j=1
|A i,j | ≤ d tr(A).

Proof. Let sgn(A) be the matrix of entries (sgn(A i,j )) 1≤i,j≤d . Then ||A|| L 1 = A : sgn(A). In addition, sgn(A) is a symmetric diagonalizable matrix and its eigenvalues α i , 1 ≤ i ≤ d, satisfy

|α i | ≤ max X∈R d sgn(A)X • X ||X|| 2 .
We have, since the entries of sgn(A) are either 1 or -1,

|sgn(A)X • X| = d i=1 X i d j=1 sgn(A i,j ) X j ≤ d i,j=1 |X i X j | ≤ 1 2 d i,j=1 (X 2 i + X 2 j ) = d ||X|| 2 .
Thus |α i | ≤ d, for 1 ≤ i ≤ d. Since A is symmetric and positive, all its eigenvalues are non-negative and the results follows from the matrix Ky Fan inequality:

A : sgn(A) ≤ Γ(A) • Γ sgn(A) ,
where, for a symmetric matrix B, Γ(B) stands for the vector obtained from the ordered set of the eigenvalues of B.

Taking the trace of Equation ( 45 Under the constraint [START_REF] Tomé | Application of the log-conformation tensor to threedimensional time-dependent free surface flows[END_REF] for the time-step, the linear System (45) for c n+1 is regular. Thus, by continuity, since tr(c n ) < b, tr(c n+1 ) < b for δt small enough. Let us suppose that there exists some values of the time step satisfying [START_REF] Tomé | Application of the log-conformation tensor to threedimensional time-dependent free surface flows[END_REF] and such that tr(c n+1 ) ≥ b; still by continuity, we may define δt 0 as the smallest time step such that tr(c n+1 ) = b, which yields:

Since the inlet condition c L is constant, the solution of this problem is stationary for t ≥ 2. The form of the source term a(x) c is chosen to mimick the production term (∇u) c + c (∇u) t involved in the transport equation of the conformation tensor. In addition, we are interested here in the situation where a(x) takes locally large values, so as to simulate zones of high elongation in a viscoelastic flows. Our understanding of viscoelastic flows is that the stiffness of this right-hand side is actually the origin of the so-called High Weissenberg Problem, i.e. the difficulty for numerical schemes to accurately capture (or, worse, not to blow up) for moderate-to-high the Weissenberg numbers. Finally, even if it is probably not essential, we chose a(x) in such a way that c(1, t) = c L which, to some extent, may be considered to represent what happens along a closed streamline in a viscoelastic flow. The chosen source term a(x) and the corresponding stationary solution c s (x) are

a(x) =        0 -1 ≤ x < -x 0 , -β sin( π x 0 x) -x 0 ≤ x ≤ x 0 , 0 x 0 < x ≤ 1. c s (x) =        c L -1 ≤ x < -x 0 , c L exp β x 0 π 1 + cos( π x 0 x) -x 0 ≤ x ≤ x 0 , c L x 0 < x ≤ 1.
where β and x 0 < 1 are given parameters. The stationary solution is obtained for t > 2. The problem may be made stiffer when decreasing x 0 and increasing β.

Let us consider a uniform partition 0 = t 0 < t 1 < • • • < t N = T of the time interval (0, T ), and let δt = t n+1 -t n for 0 ≤ n ≤ N -1 be the constant time step. Let us also consider the uniform partition -1 = x 0 < x 1 < • • • < x K = 1 of the space interval I with h = x k+1 -x k for 0 ≤ k ≤ K -1 the constant space step. The solution of (52) is approximated by piecewise constant unknowns c n k+ 1 2 , for 0 ≤ k ≤ K -1 and 0 ≤ n ≤ N , where c n k+ 1 2 is an approximation of the mean value of c(x, t n ) over (x k , x k+1 ). The scheme is initialized by c 0

k+ 1 2 = 1, for 0 ≤ k ≤ K -1.
Then the unknowns are updated in time by a second order in time operator-splitting technique, which reads, in a semi-discrete setting Advection I -Solve for c n+

In the ODE step, the real number a k+ 1 2 is defined as the mean value of a(x) over (x k , x k+1 ). This step is solved quasi-exactly by the same as for the upwind scheme (which, in fact, does not matter in the application case below, since the solution is almost constant near the beginning and end of the interval I).

(iv) The same MUSCL-type scheme applied to the logarithmic formulation of choice (ii).

We compare the resolution of Equation (52) with these different numerical schemes. The stiffness parameters are fixed to β = 50 and x 0 = 0.1, which indeed yields a rather stiff problem. The time step is fixed to δt = 1/(2 K), so as to satisfy the stability constraints of the (explicit-in-time) MUSCL schemes. Results for K = 200 are plotted in Figure 24, together with a convergence study with respect to the number of cells. In both cases, we consider the results at t = 10, when the solution is stationary both in the continuous and discrete cases. The first conclusion is that, even with the accurate solution of the ODE step implemented here, large errors may be obtained. Second, we observe that switching to the logarithmic formulation in the transport steps dramatically improves the scheme accuracy. On the opposite, the MUSCL discretization does not bring a decisive gain in accuracy. These conclusions motivate our choice for the scheme developed in this paper for the transport of the conformation tensor: we use the logarithmic formulation for the transport steps and keep an implicit linear and monotone (so first order upwind) scheme for their solution.

Note that, for one-dimensional problems as considered in this section, more sophisticated methods to cope with stiff source terms are proposed in the literature such as the so-called "well-balanced" schemes, designed to preserve at the discrete level specific "equilibrium" (i.e. steady) states [START_REF] Botchorishvili | Equilibrium schemes for scalar conservation laws with stiff sources[END_REF][START_REF] Greenberg | A Well-Balanced Scheme for the Numerical Processing of Source Terms in Hyperbolic Equations[END_REF][START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF]. Such schemes are not tested here because their extension to multi-dimensional flow problems is not straightforward.

Figure 1 .

 1 Figure 1. Primal and dual meshes for the Rannacher-Turek elements.

  For any symmetric positive definite tensor c, both c -ln(c) -I d and c + c -1 -2I d are symmetric positive definite [14, Lemma 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 τ s (u) : ∇u is non-negative. Equation (

Figure 4 .

 4 Figure 4. Boundary conditions, x-component of the velocity.

5. 1 . 1 .

 11 Oldroyd-B model at W i = 1 and W i = 3

Figure 5 .

 5 Figure 5. Evolution of velocity and tr(c) fields with the Weissenberg number. Lid-driven cavity, Oldroyd-B model. From top to bottom: W i = 1, W i = 3 and W i = 5 (with RT discretisation). Left: velocity magnitude and streamlines. Right: trace of the conformation tensor (in log-scale). Results obtained with mesh M4.

5 - 1 Figure 6 .Figure 7 . 5 - 1 Figure 8 .Figure 9 .Figure 10 . 5 - 1 Figure 11 .Figure 12 . 5 - 1 Figure 13 . 5 - 1 Figure 14 . 5 - 1 Figure 15 . 5 - 1 Figure 16 .Figure 17 .

 5167518910511112511351145115511617 Figure 6. Profiles of u x , u y and tr(c) along different lines for different meshes. Lid-driven cavity, Oldroyd-B model, W i = 1.

Figure 18 .

 18 Figure 18. Flow past a confined cylinder -Geometry outline and coarse version of the mesh. Minimum cell size of meshes of type MC Mesh MC1 MC2 MC3 MC4 MC5 δy min 8.0 × 10 -04 4.0 × 10 -04 2.0 × 10 -04 1.0 × 10 -04 5.0 × 10 -05 (a) MC-type meshes Minimum cell size of meshes of type MR Mesh MR1 MR2 MR3 MR4 MR5 MR6 MR7 δy min 6.0 × 10 -04 3.0 × 10 -04 1.5 × 10 -04 8.0 × 10 -05 4.0 × 10 -05 2.0 × 10 -05 1.0 × 10 -05 (b) MR-type meshes

Finally

  , we consider the FENE-CR fluid at W i = 0.7 for different values of the elongation parameter b.

Figure 22 Figure 19 .

 2219 Figure 19. Velocity magnitude, pressure and tr(c) (in log-scale) fields. Flow past a confined cylinder, Oldroyd-B model, W i = 0.7. Top to Bottom: velocity magnitude, pressure and trace of the conformation tensor (in log-scale).

  Profile zoomed in the wake of the cylinder

Figure 20 .

 20 Figure 20. Profile of first normal polymer stress τ xx along the curve coordinate s. Flow past a confined cylinder, Oldroyd-B model, W i = 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.

Figure 21 .

 21 Figure 21. First normal polymer stress τ xx profile zoomed in the wake of the cylinder. Flow past a confined cylinder, Oldroyd-B model, W i = 0.7. Initial meshes (left) and refined meshes (right).

Figure 22 .Figure 23 .

 2223 Figure 22. Trace of the conformation tensor normalised by the elongation parameter b (in log-scale). Flow past a confined cylinder, FENE-CR model, W i = 0.7. Top to Bottom: b = 10, b = 100 and b = 1000.

  ), we get1 δt tr(c n+1 ) -tr(c n ) = 2 ∇u : c n+1 -b λ b -tr(c n ) tr(c n+1 ) -d .

1 δt 0 b

 10 -tr(c n ) = 2 ∇u : c n+1 -b λ b -tr(c n ) (b -d).We thus have1 δt 0 b -tr(c n ) + b (b -d) λ b -tr(c n ) ≤ 2 ||∇u|| L ∞ ||c n+1 || L 1 , so, thanks to Lemma B.1, 1 δt 0 b -tr(c n ) + b (b -d) λ b -tr(c n ) ≤ 2 d ||∇u|| L ∞ tr(c n+1 ) = 2 d b ||∇u|| L ∞ .

2 for k ≥ 1 and c n+1 0 = 2 ,Figure 24 .

 20224 Figure 24. Comparison of the different methods of discretization of the advection step with β = 50 and x 0 = 0.1. Left: Stationary states obtained for K = 200. Right: Relative errors (in discrete L 2 -norm) with respect to the analytical solution as a function of the number of cells.

  1 3 :∂ t c + ∂ x c = 0 on I × (t n , t n + δt 2 ), c(x, t n ) = c n for x ∈ I. = c k+ 1 2 (t n + δt): ∂ t c k+ 1 2 (t) = a k+ 1 2 c k+ 1 2 (t), for 0 ≤ k ≤ K -1.

				(53a)
	ODE -Set c k+ 1 2 (t n ) = c	n+ 1 3 k+ 1 2	and solve for c	n+ 2 3 k+ 1 2

(53b)

Advection II -Solve for c n+1 :

∂ t c + ∂ x c = 0 on I × (t n + δt 2 , t n+1 ), c(x, t n + δt 2 ) = c n+ 2 3 for x ∈ I.

Appendix A. Positivity of a (discretely) transported tensor

The aim of this section is to prove that, with a first-order backward upwind scheme for the transport term, the solution of the transport equation for the conformation tensor (21d) remains definite positive, under a constraint for the time-step. We proceed in two steps. We first deal with the homegeneous transport (i.e. the sub-problem obtained by setting the reaction terms to zero), then with the pure reaction equation (i.e. the sub-problem obtained by setting the convection velocity to zero), and finally merge the arguments to obtain the positivity result for the original equation.

We suppose that the fraction g(c)/λ is equal to a constant, which we denote by 1/λ for short. However, the following arguments only need that g(c)/λ be a known and non-negative value in each cell of the mesh. Hence, they also hold for the FENE-CR model if the time-discretization is such that the coefficient g(c)/λ is computed with the beginning-of-step value of c, c n , provided that tr(c n ) is lower than b in any cell of the mesh. This is indeed the case with the scheme used in this paper, at the price of constraints for the time step (see Section B), which will need to be satisfied as the same time as those appearing in this section.

A.1. Pure transport

We first recall the arguments showing the positivity of the conformation tensor in the continuous setting, before mimicking the proof in the discrete setting.

Let c(x, t) ∈ R d×d be a tensor satisfying

with c(x, 0) positive and definite. The problem is supposed to be posed over a finite domain Ω ⊂ R d and, for the sake of simplicity, we suppose that u • n ∂Ω = 0 over the whole domain boundary and that the velocity field u is divergence-free. Let X ∈ R d be a (fixed) vector. Then, for 1 ≤ i, j ≤ d, we have

Summing over the indices i and j yields

Since the transport equation satifies a maximum principle and that (c(x, 0)X, X) > 0 over Ω by assumption on the initial data, we get that (c(x, t)X, X) > 0 over Ω × R, i.e. that the tensor c remains positive definite at all locations and times. The following lemma states a discrete counterpart to this property.

Lemma A.1. The discretization of (37) by the first-order backward-Euler upwind scheme preserves the positivity of the tensor c.

Proof. The discretization of ( 37) by the first-order backward-Euler upwind scheme reads

where u K,σ is an approximation of u • n K,σ such that the following divergence-free constraint is satisfied (remember that u is supposed to be divergence-free):

and

a unit eigenvector associated with this zero eigenvalue. We have, by Equation ( 42),

which raises a contradiction since, by assumption, (c n X, X) > 0.

A.3. The complete transport-reaction equation

Let us now consider the complete transport equation for the conformation tensor:

where the velocity field is supposed to be discretely divergence-free (i.e. satisfies Equations ( 38)) and the conformation tensor value at the face σ, c n+1 σ , is given by the first-order backward upwind choice [START_REF] Lozinski | An energy estimate for the Oldroyd-B model: theory and applications[END_REF]. We introduce the following set of constraints for the time step:

We are now in position to state the following result.

Lemma A.4. Let the time step be such that Equations (44) are verified, and let us assume that, for K ∈ M, c n K is a symmetric positive definite tensor. Then the linear System (43) is uniquely solvable and, for K ∈ M, the tensor c n+1 K solution of (43) is symmetric positive definite.

Proof. The proof of this result is obtained by combining the arguments of the two previous sections. First, for each cell K, the upwind choice [START_REF] Lozinski | An energy estimate for the Oldroyd-B model: theory and applications[END_REF] implies that the convection operator produces, in the system for the tensor components associated with K, an additional non-negative term on the diagonal and non-positive off-diagonal additional entries. Moreover, thanks to Relation [START_REF] Lorin | Convection systems with stiff source terms[END_REF], the diagonal additional term is equal to the sum of the off-diagonal ones. Inequality [START_REF] Renardy | Mathematical analysis of viscoelastic flows[END_REF] thus still implies that the diagonal of the matrix of this system is dominent.

Then, as in Section A.2, let us multiply the right-hand side of (43) by α ∈ [0, 1]; we easily check that the obtained system is still uniquely solvable, and denote by (c α K ) K∈M its solution. By Section A.1, we know that, for α = 0, c α K is symmetric definite positive, for any K ∈ M. If, for α = 1, c α K is not definite positive in at least one cell, by continuity, there exists α 0 ∈ (0, 1] and K 0 ∈ M such that c α0 L is non-negative in all cells L ∈ M and c α0 K0 has a zero eigenvalue. Let X ∈ R d be a unit eigenvector associated with this zero eigenvalue. We have, by the same computations as in Sections A.2 and A.1,

and, since c n K0 is definite positive, (c α0 K0 X, X) = 0 and (c α0 L X, X) ≥ 0 for L ∈ M, we easily check that the upwind choice for c α0 σ implies that the left-hand side is negative, which raises a contradiction.

Appendix B. Stability of the ODE step for the FENE-CR model

We study in this section the stability of the following semi-implicit scheme for the ODE step, for the FENE-CR model:

We observe that the two terms at the left-hand side are positive. This inequality is thus impossible, and we have tr(c n+1 ) < b for any value of the time step, if

Let us suppose that this condition is not satisfied. Inequality (47) also yields:

which, by contradiction, shows that tr(c n+1 ) < b provided that

Finally, keeping all the terms of the left-hand side of ( 47), the same argument yields that tr

The right-hand side may be seen as a function of the quantity b -tr(c n ) defined on the part of R where the denominator is positive, and which takes a minimum value. This latter thus gives an upper bound for δt which is independent of c n and reads:

Summing up the developments of this section, we get the following result.

Appendix C. A stability property of the upwind convection of tensors

Let c ∈ R d×d , u be a divergence-free regular velocity field and let F : R d×d → R be a regular function. Then, we have

The aim of this section is to state and prove a discrete analogue of this relation with an implicit upwind discretization of the convection operator. This result turns out to be an inequality and holds only when the function F is convex, as stated in Lemma C.1 below. A similar result for a scalar unknown is an easy consequence of [START_REF] Herbin | On some implicit and semi-implicit staggered schemes for the shallow water and euler equations[END_REF]Lemma A2]. The specialization of Lemma C.1 to the function of interest in the stability analysis of the scheme, namely c → -tr(ln c), is used and proved in [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF].

By convexity of the function F, we have T ∂/∂t ≥ α F(c) -F(c ⋆ ) . In addition, thanks to the equality (51), we may recast T div as

By convexity of F, R σ ≤ 0. In addition, when u σ ≥ 0, we have c σ = c (by assumption) and R σ vanishes. The product u σ R σ is thus non-negative, which concludes the proof.

Let us now apply this lemma to provide the argument requested in Section 3.2 for the stability of the scheme, i.e. a discrete analogue of Relation [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations[END_REF]. First, we recall the following classical results for functions of tensors [START_REF] Boyd | Convex Optimization[END_REF]. This is the identity needed to prove the stability of the scheme.

Appendix D. A scalar model problem

We consider the following model problem on I × (0, T ), with I = (-1, 1), ∂ t c + ∂ x c = a(x) c on I × (0, T ), c(x, 0) = 1 for x ∈ I, c(-1, t) = c L for t ∈ (0, T ). (52)