A STAGGERED PROJECTION SCHEME FOR VISCOELASTIC FLOWS

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


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 (, , ) in the following coupled system of equations with   =   (∇ + (∇)  ) and where Ω is an open bounded connected subset of R  with  ∈ {2, 3}, (0,  ) a finite time interval with  > 0,  the velocity field of the fluid,  the pressure, and  the symmetric positive-definite polymer conformation tensor.The operator  • ∇ is defined by ( • ∇)  =  • ∇  = div(  ) or ( • ∇) , =  • ∇ , = div( , ), 1 ≤ ,  ≤ , for  a generic vector or tensor field, respectively.The constant coefficients   ,   , , and  are respectively the solvent viscosity, the polymer viscosity, the fluid density and the polymer retardation time, and   stands for the identity matrix of R × .The coefficient  serves to artificially eliminate  • ∇, so as to obtain the unsteady Stokes equations for  = 0 and the Navier-Stokes equations for  = 1.The tensor   is the part of the stress accounting for the fluid solvent and   is the part of the stress accounting for the presence of polymers.The functions  () and () vary depending on the viscoelastic constitutive equations (see [10,45] for a review) with, for example, the Oldroyd-B model given by  () = () =   and the FENE-CR model by  () = () =  −tr()   , with  a real number greater than the space dimension .The dimensionless parameters that characterize the flow are the Reynolds number, Re =   /(  +   ), and the Weissenberg number, Wi = /, where  and  are a characteristic velocity and length scale.This system must be complemented by initial conditions for the velocity and the conformation tensor and by suitable boundary conditions.
A variety of numerical schemes have been proposed for the resolution of this system (see [2] for a recent review).All of these schemes face the long-standing high-Weissenberg number problem (HWNP), which makes numerical simulations extremely difficult as soon as the Weissenberg number gets close to 1.In practice, steep gradients of the stress field rapidly develop, making it difficult to correctly capture the profile of the conformation tensor and to guarantee its positivity.One of the most significant advances in pushing back the limits of the HWNP was the log-formulation proposed in [22,23] and used in many works [18,20,42,43,46,47].It consists in reformulating the problem using a logarithmic transformation of the conformation tensor.The idea described by the authors is that the conformation tensor varies exponentially and it is therefore easier to approach the logarithm of the conformation tensor with polynomial approximations than the tensor itself.Further, this formulation naturally maintains the positive-definiteness of the conformation tensor.It has later given birth to a variety of similar approaches, such as the square-root transformation [4] or the kernel conformation transformation [1].Despite these advances, the HWNP also often limits applications to moderate Weissenberg numbers because of computational costs and much progress is still needed.
One of the most important challenges of the HWNP can be understood by considering the one-dimensional transport of the scalar  posed, say, on (0,  1 ) × (0,  ),    +    = (), (, 0) =  0 (),  ∈ (0,  1 ) and (0, ) =   ,  ∈ (0,  ), where () is here to mimick the velocity gradients.This problem is similar to that of [23] but with the difference that we consider here that () can become very large and sharply varies in space.When this is the case, obtaining accurate solutions of this toy problem and, by extension, of the HWNP becomes particularly difficult.This toy problem actually belongs to a class of hyperbolic conservation laws with stiff source terms [17,38,39] that are found, for example, in modelling shallow-water flows with non-flat bottoms [3,21].A classical method, such as a centered discretization of the source term, yields non-physical solutions [9].The primary difficulty is to preserve the steady-state solutions in the discrete version of the operators.Schemes that verify this property are usually termed "well-balanced", as introduced in [11,26,33].One of the simplest ways to do this is to incorporate the discrete equilibrium states into the fluxes, which is the principle of the Upwinding Sources at Interfaces ("U.S.I.") schemes [9,11,12,34].These methods are much more accurate than classical approaches but seem next to impossible to generalise to non-scalar 2D or 3D problems.Some viscoelastic flows that admit a steady-state solution (for instance, the so-called lid-driven cavity problem, at least at moderate Wi numbers) may be schematically understood by considering the steady transport problem    = () with a field () generating, by a change of sign, successively opposite variations of .This may be understood as variations of the tensor  along closed streamlines.In Appendix D, we solve this steady-state problem with large values of () to mimic the HWNP.We show that a logarithmic transformation of the transport equation drastically improves the accuracy of the scheme, even though it does not preserve the steady-state and is therefore not well-balanced.It has also been suggested in [37,40] that a numerical scheme for viscoelastic flows should preserve, at the discrete level, a free-energy estimate derived at the continuum-level.For the system of equations (1), the time derivative of the Helmholtz free-energy is equal to the opposite of the dissipation terms plus energy sources from volume forcing terms or boundary conditions.For an initial-value problem with no sources and homogeneous Dirichlet boundary conditions, the Helmholtz free-energy must therefore decrease in time (see [7,8,31]).An ideal numerical scheme should preserve this behaviour at the discrete level and, in so doing, ensure the longtime stability of solutions.Following this idea, Barrett and Boyaval [5,6] and Boyaval et al. [14] developed numerical schemes that verify the free-energy estimate and proved the global-in-time existence of solutions to those discretizations, together with convergence to weak solutions of the continuous problem, under time-step restriction or with an additional dissipative term in the conformation equation.The authors also suggest in [13] (see also [41]) that this approach is not sufficient to guarantee accuracy.
Here, we aim for a numerical scheme that (1) preserves the positive definiteness of the tensor c, (2) satisfies a discrete equivalent of the the free-energy estimate, (3) minimizes CPU costs and is well suited to high-performance computations, (4) shows a good accuracy, in particular to capture steady-state solutions, (5) can be readily used for a variety of viscoelastic constitutive laws.
We present in Sections 2 and 3 an algorithm that satisfies Point (1) and, for the Oldroyd-B model, Point (2).We use a piecewise constant space approximation for the conformation tensor, which in fine (and especially for the final algorithm, see next paragraph) offers greats advantages from an algorithmic point of view.Equation (1c) is then solved using a finite-volume scheme.To ensure the so-called inf-sup stability while keeping a consistent order of accuracy (and so, of CPU cost), we combine this discretization with a staggered-in-space approximation of the velocity/pressure pair, with either the Marker-And-Cell (MAC) scheme [28,29], or a nonconforming low-order finite element, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedric cells (extension to simplicial cells with the Crouzeix-Raviart element [19] is straightforward).To keep the algorithm efficient and robust, we introduce a fractional-step algorithm where the solution of the Navier-Stokes equations (1a) and (1b) is computed by a standard projection method.Satisfying the discrete free-energy estimate has however a price: the transport equation for  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 , 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   .This is in fact natural for finite element discretizations, but not so for the MAC scheme for which this discretization is thus derived in Section 3.3.We also observe and prove that this weak formulation enjoys a built-in Lax-Wendroff weak consistency property [24].
Unfortunately, respecting Points (3) and ( 4) requires that we make a compromise regarding the time discretization, as detailed in Section 4. We proceed by decoupling the system of equations for the projection step, using an explicit-in-time approximation of the conformation tensor in the velocity correction equation, which greatly reduces the cost of each iteration but breaks the proof of the free-energy estimate.Consequently, Point (2) is not strictly satisfied anymore: the space semi-discretization of the momentum balance equation is consistent with the free-energy estimate, but not the time discretization.For accuracy, we further introduce a log transformation of the (now decoupled) transport equation for the conformation tensor.Unlike the standard approach of [23], however, this transformation is only partial.In the spirit of [43], we proceed by operator splitting in equation (1c), with a change of variable from  to log() 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  and, in contrast with [43], these ODEs are solved directly for , and not log(), so as to avoid any artificial introduction of nonlinearities, together with a cumbersome factorization of the velocity gradient tensor.The viscoelastic constitutive law is only involved in this cell-by-cell procedure (up to an explicit computation of   for the velocity prediction step), so that it can be easily modified to capture a variety of viscoelastic behaviors (Point ( 5)).To further gain in efficiency, we use a local time step for the solution of the ODE on each mesh element; this allows us to adjust the time step for each cell to a value small enough to preserve the positive definiteness of  (so Point (1) is still satisfied, see Appendix A.2 for the proof) while preventing a blow-up of the CPU cost.Indeed, in practice, the velocity gradient is large only in a small number of cells where the ODE becomes stiff and subcycling is needed.This greatly limits the computational cost, although one must be careful in the domain decomposition for parallel computations to keep the computing charge of the processors equilibrated.Summing up, our approach is thus based upon a low-order numerical scheme that is very efficient and reasonably accurate, the most important cause of inaccuracy being tackled by the log transformation in the transport step and a time-step adaptation to cope with the stiffness of the reaction terms.The constitutive laws for the viscoelastic fluids is essentially isolated in the resolution of cell-by-cell ODEs and can be easily modified.The resulting solver is available and open-source, as a specific module of the CALIF 3 S platform [16] developed at the French Institut de Radioprotection et de Sûreté Nucléaire (IRSN).
In the last part of the paper, Section 5, we study the efficiency and accuracy of the scheme using two benchmark 2D tests, the lid-driven cavity (up to a Weissenberg of 5) and the cylinder confined in a channel.Using massively parallel computations, we obtain solutions for the cylinder problem that are, to our knowledge, more accurate than already published results.In addition, we take benefit of the versatility of our algorithm with respect to the viscoelastic constitutive law to perform a systematic comparison of the results obtained with the Oldroyd-B and the FENE-CR models.

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: -quadrilaterals ( = 2) or hexahedra ( = 3) for a general polygonal or polyhedral domain Ω, -rectangles ( = 2) or rectangular parallelepipeds ( = 3) for a domain Ω with boundaries that are hyperplanes normal to a coordinate axis.In this case, the mesh is defined by a union of Cartesian non-uniform grids.
The sets E and E() are, respectively, the set of all ( − 1)-faces  of the mesh and of the element  ∈ 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 the cells  and  is denoted by  = |.The vector  , is the unit outward normal vector to a face  of .For 1 ≤  ≤ , we denote by ext the subset of the faces of E, E int and E ext , respectively, that are perpendicular to the th unit vector of the canonical basis of R  (these definitions being useful only in the case of Cartesian meshes).For  ∈ M and  ∈ E, || is the measure of  and || the ( − 1)-measure of the face .
The space discretization is staggered, using either the Marker-And Cell (MAC) scheme [28,29] or nonconforming low-order finite element approximations, namely the Rannacher and Turek (RT) element [44] for quadrilateral or hexahedral meshes.The degrees of freedom for the pressure and the conformation tensor (i.e. the discrete pressure and conformation tensor unknowns) are associated with the cells of the mesh M and are denoted by {  ,   ,  ∈ M}.For the velocity, we have the following cases: -For the Rannacher-Turek discretization, the discrete velocity unknowns are located at the center of the faces of the mesh and we choose the version of the element where they represent the average of the velocity through a face.For the sake of simplicity, we suppose in this presentation that the velocity is prescribed to zero over the whole boundary of Ω and these Dirichlet boundary conditions are taken into account by setting the velocity unknowns associated with an external face to zero.The set of discrete velocity unknowns reads -For the MAC discretization, the degrees of freedom for the th component of the velocity are located at the centre of the faces  ∈ E () , so, with the same boundary conditions, the whole set of discrete velocity unknowns reads { , ,  ∈ E int , 1 ≤  ≤ }.Hence, there are  unknowns per face of the primal mesh in the case of the RT scheme, corresponding to the  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  ∈ M is a rectangle or a cuboid, we define  , , for  ∈ E(), as the cone with basis  and with vertex the mass center of  (see Fig. 1).We thus obtain a partition of  in  sub-volumes, where  is the number of faces of the cell ( = 4 in two space dimensions and  = 6 in three space dimensions), each sub-volume having the same measure | , | = ||/.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  , is referred to as the half-diamond cell associated with  and .Note that the shape of  , does not need to be specified (and, for highly distorted meshes, if || is too small,  , cannot be a cone).For  ∈ E int ,  = |, we now define the diamond cell   associated with  by   =  , ∪  , .-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  ∈ M and  ∈ E(), is now the rectangle or rectangular parallelepiped of basis  and of measure | , | = ||/2 (see Fig. 2).
We denote by |  | the measure of the dual cell   and by  =   |  ′ the face separating two diamond cells   and   ′ .The set of the (dual) faces of   is denoted by Ē(  ).
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 () S associated with the degrees of freedom of the th component of the velocity (S stands for "scheme"): int for the MAC scheme, E int for the RT scheme.

A free-energy preserving pressure correction scheme for Oldroyd-B fluids
Let us first recall the formal derivation of the energy estimate given in [31].Here, we restrict the exposition to the Oldroyd-B model, and the starting point is system (1) where we set  () = () =   .We consider that no-slip conditions are prescribed on the whole boundary.Consequently, there is no source of energy in the system, either through volume forces or boundary conditions, so that we expect a decrease of the free-energy due to dissipation processes.The first step to show this consists in taking the inner product of the momentum balance equation (1a) with , which yields where  kin is the kinetic energy defined as Let us also consider the trace of equation (1c) multiplied by   /(2 ).In doing so, we suppose that any tensor  solution of equation (1c) is symmetric definite positive.We have the relationships tr(∇ ) = ∇ :  and tr( (∇)  ) = ∇ : , thanks to the fact that tr() = tr(  ) and tr( ) =  :   for any symmetric matrix  of R × and any matrix  of R × .We obtain Since  is supposed to be invertible, we may also multiply equation (1c) by  −1 to get We are now going to consider the trace of this equation.For the first term on the left-hand side, we know (see [14], Lem.1.2) that tr( −1   ) =   tr(ln()).For the second term on the left-hand side, consider that with    the tensor obtained by deriving each entry of  with respect to the th space coordinate.By the same identity as for the first term, we therefore have tr( −1  • ∇) =  • ∇tr(ln()).For the first term on the right-hand side, we remark that tr( −1 ∇ ) = tr(  −1 ∇) = tr(∇) = div() = 0 since tr( ) = tr( ) for any matrix  and  of R × .For the second term on the right-hand side, we also have tr((∇)  ) = 0. Therefore, taking the trace of equation ( 5) and multiplying by   2 , we get which can also be written as considering that We can now build a transport equation for the free-energy by adding equation (3) to equation ( 4) and subtracting equation ( 7) to obtain where , the free-energy, and   are defined by Finally, integrating over Ω yields d d by integration by parts and using the no-slip boundary condition.For any symmetric positive definite tensor , both −ln()−  and + −1 −2  are symmetric positive definite ( [14], Lem.1.1).Thus, tr(+ −1 −2  ) > 0, tr( − ln() −   ) > 0 and the previous relation shows that the free-energy  decreases with time, since the viscous dissipation   () : ∇ is non-negative.Equation ( 9) thus provides a stability estimate for the solutions of (1).

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

Prediction step -Solve for ũ𝑛+1 : Correction step -Solve for  +1 and  +1 : This time semi-discretization satisfies the stability estimate stated in the following result.
Proposition 3.1.Let (  ,   ,   ) 0≤≤ be the solution of System (10).For 0 ≤  ≤  , let us suppose that   is positive definite and let the discrete free-energy of the solution be defined by Then we have the following stability estimate Consequently, the integral over the domain of the discrete free-energy is non-increasing.
Proof.Using the identity 2  ( − ) =  2 + ( − ) 2 −  2 valid for any pair of real numbers  and , we have ũ+1 Hence, using the fact that   is divergence-free to recast the transport term as a convection term, thanks to the no-slip boundary condition, Thanks once again to the boundary conditions, we have Therefore, taking the inner product of equation (10a) with ũ+1 and integrating over the computational domain yields Let us now rewrite equation (10b) as Taking the square of the  2 -norm of both sides and multiplying by /(2 ) yields 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) Adding with (15), we obtain the following discrete kinetic energy inequality with, for  =  and  =  + 1, Let us now multiply (10d) by For any symmetric definite positive matrices  and , we have tr( −1 ( − )) ≤ tr ln() − tr ln() (see [14], Lem.2.1), which yields Taking the trace of relation (18) and integrating over Ω, by arguments already used in the continuous setting, we thus get Taking the trace of equation (10d), adding the discrete kinetic energy balance (17) and subtracting the latter relation as in the continuous setting yields the desired estimate, thanks to the duality relation

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 ũ+1 : Correction step -Solve for  +1 ,  +1 and  +1 : In the prediction step and in the velocity correction step, for  =  and  =  + 1, the tensor    is discretized on the primal cells and is computed as a function of the conformation tensor as 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  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 ∇ +  (∇)  instead of ∇  +  (∇)  ), as is done in equation (21d).Once the equation is symmetrized, if   is symmetric, for 1 ≤ ,  ≤ , the discrete equation for  +1 , is strictly the same as the equation for  +1 , .Conversely, a symmetric tensor solution to (21d) is also a solution of the original equation and thus switching to (21d) preserves the consistency of the scheme.Moreover, it is proved in Appendix A that a first-order backward upwind scheme for (21d) preserves the positive definiteness of , under a condition for the time-step depending on the  ∞ -norm of ∇ +1 .Finally, for this equation and this convection scheme, we need a discrete analogue of (19) which is proved in Appendix C.
Let us now turn to the momentum balance equation.A finite-volume discretization of the convection term for the RT discretization may be found, in the more general context of variable density flows, in Section 4.2 of [35].The discrete kinetic energy balance (i.e. the discrete equivalent of ( 13)) is established in Lemma 7.2 of [35].For the MAC scheme, a discretization of the convection term is proposed in Section 2, equation ( 29) of [24] and its  2 -stability (i.e. as before, the discrete equivalent of ( 13)) is shown in Lemma 3.6 of [24].The convection operator takes the following general form where  , is an approximation of the mass flux through  computed from the fluxes through the primal faces (i.e. for  ∈ M and  ∈ E(), the quantities  , =  ||  , with  , given by (24a) and ( 29) for the RT and MAC schemes respectively) in such a way that a divergence-free constraint is also satisfied on the dual cells   (i.e.∑︀ ∈ Ē()  , = 0).A construction of the dual mass fluxes may be found in Appendix A of [36] (RT approximation) and Section 2 of [24] (MAC scheme).The approximation of the th component of the velocity at the dual face, (ũ +1  )  , may be chosen upwind or centered.Finally, we need to ensure the duality of the velocity divergence and pressure gradient, together with the duality relations ( 14), ( 16) and (20).All these properties are obtained "by construction" when the associated terms in the momentum balance equation are derived from a variational formulation of this latter equation, as is done with the RT discretization (see Sect. 3.3.1 below).For the MAC scheme, the pressure gradient-velocity divergence duality is standard.The definition of these operators may be found, for instance, in Section 2 of [24] and a proof of their duality is given by Lemma 2.4 of [24].For the viscous diffusion term, we propose in [25] a construction that ensures the positivity of the dissipation, i.e. a discrete analogue of Relation (14).It consists in obtaining the usual finite-volume discretization of this term by a variational technique, thanks to a discrete  1 -inner product defined for "finite-volume discrete functions".For the sake of completeness, we recall these developments in Section 3.3.2below, 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 ,   , 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 , ( ) of the total stress tensor  = −  +  +  .We want this quantity to satisfy a discrete counterpart of the duality relations ( 14), ( 16) and ( 20) that, as shown in the previous sections, are crucial for the preservation of the free-energy estimate by the scheme.Gathering these three relations in a single one, we thus search for an identity of the form: where  M (, ) stands for non-negative mesh-dependent bilinear form associated to the viscous dissipation and the terms div   and ∇  , to be used in the discretization of equations ( 21c) and (21d) respectively, have to be defined.The derivation of div ,  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 d, where stands for the vector-valued finite element shape function associated with the th component of the velocity and the face .By definition of the RT finite elements, this shape function is    () , 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  1 -mapping from the space spanned by the set of functions } on the reference (0, 1) 3 cube ( = 3).By a simple reordering of the sums, equation (22) reads Considering separately the pressure gradient, the diffusion and the elastic extra-stress terms in  , we obtain the following expression for the discrete operators featured in the scheme (21), dropping the time index for short with For 1 (24e)

MAC scheme
For the MAC scheme, the strategy implemented here is to mimick the RT case, i.e. to recast the momentum balance equation under a weak form.We define a discrete gradient operator acting on the velocity discrete functions (i.e.piecewise-constant functions over the dual cells).Each component of this discrete gradient is itself a piecewise-constant function, based on a new partition of the computational domain.In two-dimensions, the discrete analogues of     and     are defined over the primal cells, while the other components of the velocity gradients are defined over volumes associated with the mesh vertices.With this gradient operator defined, the machinery of standard Galerkin techniques applies for the discretization of the momentum balance.The discrete velocity function space is shown to be spanned by a set of shape functions, and each discrete equation is derived by using these shape functions as test functions.Thanks to this construction, the scheme discrete operators natively satisfy the duality Relation (22) (they are indeed obtained from this relation itself), and have formally the same expression than in the finite element case, i.e. obey formula similar to System (24) above.
This approach was followed in [25] for Newtonian fluids.We briefly recall here these results and use the same technique to cope with the divergence of the elastic extra-stress tensor.For short, we only address here the two-dimensional case and we refer to [25] for the (natural) extension to the three-dimensional case.
The discrete velocity gradient.The developments of this section rely on the geometric properties of structured meshes and we need to momentarily switch to the notations that are specific to this context.We thus suppose that the primal mesh reads M = { , , (, ) ∈ I}, with I ⊂ N 2 with notations for geometrical quantities that are given in Figures 3-5.We detail here the discretization of the terms associated with the first component of the velocity, the adaptation to the second component being straightforward.In the interior of the computational domain, the discrete partial derivatives of this velocity component are defined as follows.
-For (, ) ∈ I, let the primal cells be ).The discrete derivative involved in the divergence (so, for the velocity -component, only  M    ) is defined over the primal cell by, ∀(, ) ∈ I and ∀ ∈  , , -For the other derivatives (so, for the velocity -component, only  M    ), we introduce a fourth mesh which is vertex-centred, and we denote by   the generic cell of this new mesh, with   ( −1 ,   ).We define I  the set of pairs (, ) such that  − 1 2 ,− 1 2 is an internal vertex of the mesh (i.e. a vertex which does not lie on the boundary).Then, for (, ) ∈ I  and ∀ ∈   − 1 2 ,− 1 2 : The only necessary extension of this definition to cope with boundaries concerns the definition of  M    over a half vertex-centered cell associated with a vertex lying on a horizontal boundary (see Fig. 5).In this case, we use the usual "fictitious cell trick" in order to apply Relation ( 26): an external cell, of zero -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 -component of the velocity, the discrete diffusion tensor can be defined as Finite-volume test functions.Let us denote by ) is the mass center of a vertical (resp.horizontal) internal face of the mesh.For (, ) ∈ I  , we denote by  ,(− 1 2 ,) the test function associated with the degree of freedom of the -component of the velocity located at  − 1 2 , .This discrete function is defined as Its non-zero partial derivatives are  M   ,(− 1 2 ,) and  M   ,(− 1 2 ,) with , 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  with its associated piecewise constant function, we have Using the expression of the discrete derivatives of the test functions, we get 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 which is a definition similar to the RT discretization, equation (24b).The velocity divergence operator immediately follows from the definition of the velocity discrete partial derivatives so, since the partial derivatives of the velocity involved in the divergence are constant over each primal cell, 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 ||  , , with, for  ∈ E () ,  , =    , •  () . ( Finally, the discrete divergence of the viscous stress tensor is defined by the following weak formulation By a standard reordering of the sums, we get An important consequence of this relation, with the definition of  M  (ũ) =   (∇ M ũ + (∇ M ũ)  ), is that, thanks to an elementary property of the tensors contraction, we may replace the second tensor ∇ M ũ in (31) by its symmetrization (∇ M ũ + (∇ M ũ)  )/2.We thus obtain that the dissipation (i.e. the integral in (31)) is nonnegative, which is a discrete counterpart of (14).In addition, the duality Relation (20) requires the discretization of ∇   in the right-hand side of the transport equation for the conformation tensor (21d) to be This weak formulation of the diffusion term (31) is shown in [25] to lead to the usual diffusion operator when the viscosity does not depend on space.The same is true (i.e. the weak formulation yields the usual conservative finite-volume scheme) for a space-dependent viscosity only when this quantity is considered as a piecewiseconstant function over the same control volumes as the discrete velocity gradient components.For instance, for two-dimensional spaces, in the compression terms ( M  ) , and ( M  ) , , the viscosity must be associated with the primal cells, while, for the shear components ( M  ) , and ( M  ) , , 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   to a piecewise function over the primary cells by Then we set ∀(, ) ∈ I  , −(div  ) 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   − 1 2 , of a discretization of the flux associated with div  , i.e. the integral of the first component of    , (see Fig. 3).
-On the vertical face  =   , , the flux is given by ∫︁ -The flux through the horizontal half face )︁ , with 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 , whereas, with a linear interpolation, we would have , 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.Lemma 3.2.Let (M () ) ∈N be sequence of structured discretizations of Ω, with a space step tending to zero.Let ( () ) ∈N be a sequence of discrete conformation tensors, with  () associated with M () .We suppose that ( () ) ∈N weakly converges in  1 (Ω) × to .Let  ∈  ∞  (Ω)  , and let us denote by  () the interpolate of  in the space of velocity discrete functions defined by  () () =   (  ) for  ∈   ,  ∈ (E () ) () , 1 ≤  ≤ , with   the mass center of .Finally, for a discrete tensor  , i.e. a tensor-valued piecewise constant function over the primal cells of M () , we denote by div M ()  : Ω → R  the vector-valued function defined by Proof.By construction, we have and the results follows thanks to the convergence of ∇ M  () () to ∇() for almost every  ∈ Ω.
This result extends to nonlinear laws, under stronger convergence and boundedness assumptions (and/or regularization of the function ) for the sequence of conformation tensors (in fine, what is needed is the weak convergence in  1 (Ω) × of ( () ) ( () −   )).

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 ,  +1  − div ,    ) tightly couples the equations of this step and thus makes its resolution time-consuming.In addition, it does not seem to bring any significant gain in accuracy.This term is thus suppressed, so as to decouple the transport equation of the conformation tensor (21d) from equations (21b) and (21c) (which, as usual, are recast under the form of an elliptic problem for the pressure).
-Second, especially at high Weissenberg numbers, the production terms in equation (21d) (i.e. the terms ) are likely to become very stiff, so that stability and accuracy of usual schemes may require prohibitively small space and time steps.The analysis of a toy scalar problem (see Sect.D) shows that the scheme accuracy is dramatically enhanced by making a change of variable from  to log().However, for general equations of state (in particular, as in the FENE-CR model, when the function () does not boil down to identity), this change of variable brings additional nonlinearities.A solution for this problem is given by the observation that, still on the toy problem, the scheme accuracy is essentially preserved by a 3-steps Strang-type decoupling, that consists in performing one half-step of homogeneous transport of log(), then dealing with the reaction terms with the original variable  and finishing by the second half-step of transport of log().Up to the choice of the variable in the middle step, a similar strategy was already adopted in [43], with a different space discretization.With the specific space discretization considered here, i.e. with piecewise-constant discrete conformation tensors, the middle step consists in a family of Ordinary Differential equations, one per primal cell, for all the tensor components.Since the stiffness of this step is directly correlated to the magnitude of the velocity gradient in the cell, this allows us to adopt a local time-stepping strategy with a time step per cell, which dramatically enhances the accuracy and robustness of the scheme.Moreover, in most applications, we observe that large velocity gradients are localized in very specific zones of the flow so that small local time steps are used only in a few cells and the scheme efficiency is only slightly affected.
Finally, we obtain a two-steps scheme, where we successively deal first with hydrodynamics and then with the conformation tensor transport.The hydrodynamics step reads Prediction step -Solve for ũ+1 : Correction step -Solve for  +1 and  +1 : Up to the (explicit) term div ,    in the prediction step, this system is the usual so-called incremental pressurecorrection 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  + 1 3 : and solve for  + 2 3 =   (  + ): Advection II -Solve for  +1 : Transport steps are discretized by a standard backward first-order upwind scheme.From the form of these equations, these steps preserve the symmetry of the conformation tensor: if   is symmetric, the same equation is solved for  +1 , and  +1 , , 1 ≤ ,  ≤ .In addition, since we search for log() and then compute  by taking the exponential of the result, these transport steps necessarily yields positive definite conformation tensors.The local ODE (33b) is solved using a first-order semi-implicit Euler scheme, with a local sub-time step   = /  , that reads, for  ∈ M and 0 with  0  =

𝐾
=     .We prove in Section A.2 that this relation yields a symmetric positive definite tensor under the constraint for the choice of the sub-time step   stated by Inequality (A.5) of Appendix A.2, which imposes a scaling of the time step as 1/||∇ +1  || ∞ .In addition, for the FENE-CR model, we prove in Section B that, under additional constraints on the sub-time step (see Lem. B.2), the solution of (34) satisfies tr( +1  ) < , for  ∈ M. Note that the constitutive relation only appears in the step (33b) and in the expression of the polymeric stress tensor in the prediction step (32a), with explicit-in-time values in the latter case; in addition, it associates a piecewise constant conformation tensor to a piecewise constant (on the same mesh) stress tensor.These two features make changes of constitutive relations very easy in practice.The ODE step (34) is already written under its general form and, for the velocity prediction step (32a), the tensor    reads: and the discretization of its divergence is described in Section 3.3.1 (RT discretization) or Section 3.3.2(MAC discretization).

Numerical results
Here we consider the dimensionless version of the system of equation (1) in the vanishing inertia limit ( = 0), where space is non-dimensionalized with  and time with / , where  and  stand for the characteristic length scale and velocity respectively, Wi = /,  =   /(  +   ) is the viscosity ratio and Re is the Reynolds number defined by Re =   /(  +   ).The pressure is non-dimensionalized with a reference pressure equal to (  +   ) /.The constitutive relation is given by ,  > , for FENE-CR.
The momentum balance equation (35a) may be complemented by a body force; in this case, it is explicitely mentioned.In all this section, we solve equation (35) with Re = 1 (remember that Re multiplies the velocity time derivative term only, the momentum convection term being removed in ( 35)), so the data of , Wi, the geometry and the initial and boundary conditions completely defines the problem.We address two benchmark cases from the literature.The first one is the lid-driven cavity, which we solve with the MAC scheme and with the RT discretization, both on Cartesian (uniform and non-uniform) grids.The second one is the flow around a cylinder confined in a channel, which we solve on an unstructured bodyfitting mesh using the RT discretization.Even though classic [18,20,22,23,27,32,43,48], each of these cases is challenging to solve numerically even at moderate Wi numbers.In the case of the lid-driven cavity, it is only with the recent introduction of the log-conformation formulation that results at Wi ≥ 1 have been obtained [22,23].For the flow around the cylinder, a singularity in the wake of the cylinder appears at Wi ≃ 0.7 that is very difficult to capture accurately [20,32].
Our goals are here to validate the proposed scheme by comparison with results from the literature, and to study its accuracy and numerical convergence properties.We present results for both the Oldroyd-B and the FENE-CR models, these results being new in the latter case, and compare the obtained solutions.In the lid-driven cavity case, we quantitatively assess the convergence of the solutions of the FENE-CR model to the solutions of the Oldroyd-B model when the parameter  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 [16].The whole domain is decomposed into several sub-domains with the METIS (4.0.3) graph partitioner and the communication between the sub-domains is done through the OpenMPI (3.1.5)interface.The linear systems are solved by Krylov methods, using a Block-Jacobi preconditioner at the global level and an ILU preconditioner for the blocks associated to the subdomains.Parallel computations were then carried out on TotalEnergies' group supercomputer Pangea II.
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 M,  = 1, . . ., 4, use a non-uniform step in the two directions.These are built from lists for the -coordinates and -coordinates of the edges, which are defined as follows.
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  = 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  = 10 −3 , which is the time step beyond which unstationary phenomena appear.The sub-time-step for the solution of the ODE on the cell  is set to /  with   the smallest integer number such that /  ≤ 1/(2  ||∇   +1 || ∞ ) for the Oldroyd-B model and /  ≤ ( − 2)/(4    ||∇   +1 || 2 ∞ ) for the FENE-CR model, with  = 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 6 the 2D representations of the velocity magnitude and the trace of the conformation tensor for Wi = 1 and Wi = 3 when the steady state is reached.The velocity field exhibits some of the classical features that have been observed in the literature.Each corner at the bottom of the cavity contains a vortex at very low velocity that effectively behaves as a stagnant zone.As we increase the Wi, we transition from a left-right symmetry (Wi = 0) to an asymmetric profile with a shift of the primary central vortex and of the region of maximum velocity close to the lid towards the left-hand side.The most interesting features are observed for the trace of the conformations tensor.We see a dramatic increase of the maximum value of the trace from Wi = 1 to Wi = 3, going from approximately 5 × 10 3 to more than 10 5 .This maximum is reached at the lid with zones of large stresses both in the middle of the lid and in the top right corner.This suggests that mesh convergence in these regions may be particularly difficult to obtain.We also observe that regions of high stresses have extended further inside the cavity at Wi = 3 compared to Wi = 1, which may impact convergence even away from the lid.
We then study mesh convergence and compare quantitatively our scheme to results from the literature, starting with the Weissenberg number equal to 1 and comparing to results in [23].The computation is performed up to a time  = 10, with a steady state reached for  ≈ 8.The first component of the velocity,   , along the line  = 0.5, is presented in Figure 7a.The steady state is almost independent from the mesh and is in close agreement with results in [23].The second component of the velocity,   , along the line  = 0.75 is plotted in Figure 7b.Convergence is slower and seems to be roughly obtained for the 256 × 256 grid, which is also in close agreement with [23].
The most difficult aspect of these computations is to accurately estimate the conformation tensor near the lid.We investigate this issue here using the M meshes that are specifically designed for this purpose, with finer meshes in the region close to the lid. Figure 7c shows tr() along the line  = 0.5 for M1, M2, M3 and M4.On this plot, convergence seems to be achieved for any of the grids.However, we can see in Figure 7d that the profile of tr() becomes very steep at the top.Figures 7e and 7f show tr() along the lines  = 0.975 and  = 1 respectively, for the different meshes.At  = 0.975, convergence seems almost achieved.The picture is completely different at  = 1 where the maximum value of tr(), obtained close to  = 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 ( + 1)/ for  = 1, . . ., 3.
We now focus on the case Wi = 3 and compare with results in [18,23].When the mesh is coarse, an unsteady behavior is observed that disappears when refining the mesh.The temporal evolution of the kinetic energy, is presented in Figure 8.For the uniform 64 × 64 grid, we observe periodic oscillations of the kinetic energy that are related to the appearance of vortices in the upper-right corner of the cavity.However, this phenomenon seems to be purely numerical, since, with finer meshes, oscillations disappear and a stationary solution is obtained after  ≈ 40.Such numerical instabilities had also been observed in [23], although they did not obtain a stable solution by refining the mesh.In [18], solutions for this case may also be slightly unsteady.We compare our stationary solution to the solutions obtained in these two papers at a given time  = 40.Figures 9a and 9b depict profiles of the velocity field for the M meshes, showing that convergence is slower than for Wi = 1.We find that our results for the first velocity component (Fig. 9a) are in very good agreement with those of   the literature [18,23].For the second component, we observe differences with both [18,23] in the region where  0.5.However, both results in [18,23] are also significantly different from one another and our solution lies roughly in between the results of these two papers.The zone where this discrepancy is observed is also the area which is likely to be impacted by unstationarities generated near the upper right corner (this was observed in our computations with coarse meshes).For the conformation tensor, difficulties in converging are no longer confined to the very top of the cavity but rather propagate further inwards (see Figs. 9c-9e).As for Wi = 1, we still obtain an acceptable convergence everywhere except in the region closest to the lid (Fig. 9f).Since momentum and conformation transport are strongly coupled, the absence of convergence for the conformation tensor locally in the lid region, which we believe is a common feature of all numerical results so far, raises concern for the convergence of the velocity field.The first clue that we may have convergence for the velocity field is of course the observed mesh convergence for the velocity components, together with the 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 Wi = 5
Finally, we consider the case Wi = 5.This case has been studied in [20,23], where the authors obtained an unsteady regime.In [23], where they use a uniform grid MAC discretization, they obtained oscillations of the kinetic energy with a frequency that decreases with the mesh size.We obtain a similar result in Figure 10 (top) with our MAC scheme on coarse meshes and a stabilization for the M1 mesh (Fig. 10, bottom).These oscillations arise from the fact that the zones where tr() is high form a very thin filament for Wi = 5 that is very difficult to capture numerically.When the mesh is not sufficiently refined, the filament can rupture, generating energy transfer and kinetic energy oscillations.Refining the mesh delays this rupture and decreases the frequency of the kinetic energy oscillations.In contrast with what happens at lower Weissenberg number, these oscillations generate slight perturbations which remain apparent on the profile of the -component of the velocity along the line  = 0.75, in the right part of the curve (let us say for 0.7 ≤  ≤ 0.9); except for the conformation tensor very close to the lid, all the other profiles do not vary with time and match the results of the RT discretization (Fig. 12, see comments below).By switching to the RT discretization, oscillations disappear and we visually obtain a plateau of kinetic energy at  ≈ 50, Figure 11.However, it seems that oscillations of very small amplitude either subsist in the numerical solution or may redevelop if the time step is relaxed.Whatever the reason for this, this prevents increasing the time step at the end of the computation.For the M4 mesh, the minimum time step is set at 10 −6 and the time step remains close to 10 −5 at the end of the computing time.However, in this case also, we believe that the mechanism leading to temporal oscillations is once again numerical.
The results obtained at Wi = 5 at steady state are shown on Figure 12 for the RT discretization.As in the previous cases, convergence is obtained on the velocity and the conformation tensor profiles, except close to the lid, for the conformation tensor.

FENE-CR model
We now study the solutions of the FENE-CR model.As we will see, convergence is more easily reached than with the Oldroyd-B model, so the MAC and RT discretizations yield equivalent results.The solutions plotted here are obtained with the MAC scheme.
The 2D plots of the trace of the conformation tensor (log-scale) for Wi = 1 and Wi = 3 are compared in Figure 13 for different values of .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 , the regions of large stresses are much more diffuse and their size tends to reduce as  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  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() is lower than .Small values of  can considerably limit the polymer elongation, the polymeric stress and the impact of polymers on flow.This is also clear when considering the limits of the FENE-CR model, with, at least formally,  →  leading to Stokes flow and  → ∞ to the Oldroyd-B model.
Figures 14 and 15 show the mesh convergence for  = 1000 at Wi = 1 and Wi = 3, respectively.For the components of the velocity, convergence is similar to that of the Oldroyd-B model.For the conformation tensor, the behavior is radically different, especially close to the lid.Convergence is much faster in the case of the FENE-CR model, as expected from the effects of the limitation of the polymer elongation.Figures 16 and 17 show results of the FENE-CR model for different values of  at Wi = 1 and Wi = 3, respectively.These are further compared with the Stokes and Oldroyd-B models.As expected, results for the FENE-CR model get closer to the Newtonian solution when  decreases and to the Oldroyd-B solution in the limit  → +∞. Figure 18 shows the evolution, as a function of , of the difference (in L2-norm) between the velocity fields: with   the velocity field obtained with the Oldroyd-B model and   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 , we converge towards the Oldroyd-B model, with a faster convergence at Wi = 1.Since results for the FENE-CR model are converged for both the velocity and conformation fields, this convergence supports the claim that the velocity profiles for the Oldroyd-B fluid are correct.

Flow past a confined cylinder
We now consider the case of a two-dimensional flow around a confined cylinder.Schematics of the geometry are presented in Figure 19.It consists of a cylinder of radius  = 1 confined in a channel of width  = 4 and length  = 30.No-slip and no-penetration boundary conditions,   =   = 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  =    (where   stands for the first vector of the canonical basis of R 2 ), with  adjusted so that the flow rate is constant and set to  = 2.We further fix  = 0.59 and perform computations for Wi = 0.6 and Wi = 0.7.
The presence of a stagnation point in the wake of the cylinder generates a thin zone of very large polymer elongation, a "filament", located on the symmetry axis and starting from the intersection of the latter with the boundary of the cylinder.This filament is often referred to as a birefringent strand in the literature, in reference to an experimental method to visualize it.It is notoriously difficult to capture in flow simulations because of both the large values of the conformation tensor and the strong gradients.Mesh convergence for the conformation tensor  and, equivalently, for the polymer stress tensor, requires a very refined discretization in this area.To solve this issue, we use a set of non-structured meshes (Fig. 19) with the RT space approximation.These meshes are built with the GMSH mesher (version 4.4) as follows: -outside of the blue zone in Figure 19 (delimited by the vertices (12.5, 0), (12.5, 0.5), (14,2), (16,2), (17.5, 0.5) and (17.5, 0), if the bottom left corner of the channel is associated with the origin), the length of the cell edges is roughly kept constant at value  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  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  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 on the cell  is set to /         Table 1.Minimum width  min of the cells located in the wake of the cylinder.

Oldroyd-B model
We present the results for two Weissenberg numbers, Wi = 0.6 and Wi = 0.7.In both cases, the overall characteristics of the flow are the same so, for two-dimensional plots, we restrict ourselves to the two-dimensional fields (velocity, pressure and trace of the conformation tensor) on the MC4 mesh for the highest Weissenberg number, once the flow has reached the stationary state (Fig. 20).For Wi = 0.7, this occurs at time  ≈ 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.To go further, we compare our results with those proposed in [20,32].For Wi = 0.6, the stationary regime is obtained at time  ≈ 2. Figure 21 shows the -component of the polymer stress, as a function of the curvilinear coordinate  along the line composed of the cylinder boundary and the part of the symmetry axis located in the wake of the cylinder.The two stagnation points at the cylinder boundary correspond to  = 0 for the front one and  =  for the rear one.As expected from the two-dimensional plot, we observe two peaks for   .The first one, more pronounced, is approximately centered at the throat between the cylinder and the solid surface,  ≈ /2, and the second one is located in the wake close to the cylinder.These results are in very good agreement with both [20,32].We now turn to the case Wi = 0.7, for which a lack of convergence on   in the wake of the cylinder is generally reported in the literature (see [20,32]).This is observed in Figure 22 (left) with the MC-type meshes.We find a profile close to the one proposed in [32] with the MC1 mesh and close to the one proposed in [20] 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 22 (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.

FENE-CR model
Finally, we consider the FENE-CR fluid at Wi = 0.7 for different values of the elongation parameter . Figure 23 shows the trace of the conformation tensor in log-scale.As for the lid-driven cavity, when the value of  increases, the zones of high elongation become thinner and profiles look qualitatively similar to those obtained with the Oldroyd-B model.
Results obtained for the -component of the polymer stress, 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 ()/ is equal to a constant, which we denote by 1/ for short.However, the following arguments only need that ()/ 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 ()/ is computed with the beginning-of-step value of ,   , provided that tr(  ) is lower than  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 Sect.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 (, ) ∈ R × be a tensor satisfying with (, 0) positive and definite.The problem is supposed to be posed over a finite domain Ω ⊂ R  and, for the sake of simplicity, we suppose that  •  Ω = 0 over the whole domain boundary and that the velocity field  is divergence-free.Let  ∈ R  be a (fixed) vector.Then, for 1 ≤ ,  ≤ , we have Summing over the indices  and  yields Since the transport equation satisfies a maximum principle and that ((, 0), ) > 0 over Ω by assumption on the initial data, we get that ((, ), ) > 0 over Ω × R, i.e. that the tensor  remains positive definite at all locations and times.The following lemma states a discrete counterpart to this property.
Lemma A.1.The discretization of (A.1) by the first-order backward-Euler upwind scheme preserves the positivity of the tensor .
Proof.Let sgn() be the matrix of entries (sgn( , )) 1≤,≤ .Then ||||  1 =  : sgn().In addition, sgn() is a symmetric diagonalizable matrix and its eigenvalues   , 1 ≤  ≤ , satisfy We have, since the entries of sgn() are either 1 or −1, Thus |  | ≤ , for 1 ≤  ≤ .Since  is symmetric and positive, all its eigenvalues are non-negative and the results follows from the matrix Ky Fan inequality: where, for a symmetric matrix , Γ() stands for the vector obtained from the ordered set of the eigenvalues of .
Taking the trace of equation (B.9), we get Under the constraint (B.10) for the time-step, the linear System (B.9) for  +1 is regular.Thus, by continuity, since tr(  ) < , tr( +1 ) <  for  small enough.Let us suppose that there exists some values of the time step satisfying (B.10) and such that tr( +1 ) ≥ ; still by continuity, we may define  0 as the smallest time step such that tr( +1 ) = , which yields: We observe that the two terms at the left-hand side are positive.This inequality is thus impossible, and we have tr By convexity of F,   ≤ 0. In addition, when   ≥ 0, we have   =  (by assumption) and   vanishes.The product     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 (19).First, we recall the following classical results for functions of tensors [15].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  × (0,  ), with  = (−1, 1),    +    = ()  on  × (0,  ), (, 0) = 1 for  ∈ , (−1, ) =   for  ∈ (0,  ).(D.16) Since the inlet condition   is constant, the solution of this problem is stationary for  ≥ 2. The form of the source term ()  is chosen to mimick the production term (∇)  +  (∇)  involved in the transport equation of the conformation tensor.In addition, we are interested here in the situation where () 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 Weissenberg numbers.Finally, even if it is probably not essential, we chose () in such a way that  Then the unknowns are updated in time by an operator-splitting technique, which reads, in a semi-discrete setting Advection I -Solve for  + In the ODE step, the real number  + 1 2 is defined as the mean value of () over (  ,  +1 ).This step is solved quasi-exactly by  otherwise.The definition of the face approximations for  = 0,  = 1 and  =  is 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 ).(iv) The same MUSCL-type scheme applied to the logarithmic formulation of choice (ii).
We compare the resolution of equation (D.16) with these different numerical schemes.The stiffness parameters are fixed to  = 50 and  = 0.1, which indeed yields a rather stiff problem.The time step is fixed to  = 1/(2 ), so as to satisfy the stability constraints of the (explicit-in-time) MUSCL schemes.Results for  = 200 are plotted in Figure D.1, together with a convergence study with respect to the number of cells.In both cases, we consider the results at  = 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 [11,26,33].Such schemes are not tested here because their extension to multi-dimensional flow problems is not straightforward.

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

Figure 2 .
Figure 2. Primal and dual cells for the MAC discretization in two space dimensions.(a) Mesh associated to the first component of te velocity.(b) Mesh associated to the second component of the velocity.

Figure 3 .
Figure 3. Unknown and dual cell for the -component of the velocity, notations for staggered discretizations.

Figure 4 .
Figure 4. Discrete partial derivatives of the -component of the velocity.

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

Figure 11 .
Figure 11.Time evolution of the kinetic energy  kin for the M meshes with RT discretization.Lid-driven cavity, Oldroyd-B model, Wi = 5.

Figure 18 .
Figure 18.Difference from the Oldroyd-B velocity field for different values of the elongation parameter .Lid-driven cavity, FENE-CR model, Wi = 1 and Wi = 3. Results obtained with the mesh M4.

Figure 19 .
Figure 19.Flow past a confined cylinder -Geometry outline and coarse version of the mesh.

Figure 20 .
Figure 20.Velocity magnitude, pressure and tr() (in log-scale) fields.Flow past a confined cylinder, Oldroyd-B model, Wi = 0.7.Top to bottom: velocity magnitude, pressure and trace of the conformation tensor (in log-scale).

Figure 21 .
Figure 21.Profile of first normal polymer stress   along the curve coordinate .Flow past a confined cylinder, Oldroyd-B model, Wi = 0.6.Left: profile along the cylinder surface and in the wake of the cylinder.Right: the same profile zoomed in the wake of the cylinder.(a) Profile along the curve coordinate.(b) Profile zoomed in the wake of the cylinder.

Figure 22 .
Figure 22.First normal polymer stress   profile zoomed in the wake of the cylinder.Flow past a confined cylinder, Oldroyd-B model, Wi = 0.7.Initial meshes (left) and refined meshes (right).(a) MC-type meshes.(b) MR-type meshes.

2 for=Figure D. 1 .
Figure D.1.Comparison of the different methods of discretization of the advection step with  = 50 and  = 0.1.Left: stationary states obtained for  = 200.Right: relative errors (in discrete  2 -norm) with respect to the analytical solution as a function of the number of cells.