D Sandri

Finite element method for viscoelastic fluids flows : approximation of the upper convected Maxwell model and study of a wall slip -Mesh convergence study

A wall slip (or rather, quasi-slip) appears in the simulation of viscoelastic fluids obeying the upper convected Maxwell (UCM) model, in the case of a flow through a 4 : 1 contraction with zero velocity conditions at the walls (cf. [Eur. J. Mech. B Fluids, 24, (2005) 733-750]). This phenomenon appears with the use of the new finite element method (FEM) that we have developed in a previous work (see [Comput. Methods Appl. Mech. Engrg. 191 (2002) 5045-5065]), the so-called θ-MSUPG method. This FEM use a splitting of the equations with a parameter θ ∈ [0, 1] and a modified upwind term in the SUPG (Streamline Upwind Petrov Galerkin) formulation. One goal of this paper is to study the influence of the splitting and of the modified upwind term in the occurrence of the slip. To do this, we perform a comparative study that requires the use of FE satisfying inf-sup Tensor-Velocity conditions. To this purpose, we use a FE inspired from a paper of V. Ruas (see [Japan Journal of Industrial and Applied Mathematics, 11, (1994), 113-130]). It appears numerically that, for the FEMs studied, the modified upwind term is necessary to obtain a wall slip phenomenon and that the use of splitting can effectively stabilize this phenomenon (see Fig. 3.19 & Tab. 5). This could explain why this wall slip is not mentioned in the literature review. A second objective of this paper is to study the mesh convergence of the θ-MSUPG method. We obtain numerically the convergence of the θ-MSUPG method towards a solution with wall slip (see Tab. 3) and we give a heuristic explanation of the values of the extra-stress tensor found at the downstream wall of the flow (see Eqn. (3.10) & Tab. 4). With mesh refinement, we observe an increase in the convergence rate of the fixed point used to solve the nonlinear FEM problem (see coefficient γ, in Tab. 3), but we did not observe the high Weissenberg number problem (divergence of the fixed point, with mesh refinment) as in the use of SUPG formulation. A recent work shows the existence of weak theoretical solutions with slip for the UCM fluid flow problem through a channel and provides an element of confirmation that the wall slip observed in the calculations of this paper is not a simulation artifact (see HAL : hal.science/hal-04088918v2). Other results are also described with simulations in a channel as well as simulations for the corotational Maxwell model (COM model) : this model gives wall slip for all types of FE methods, as soon as the imposed inflow is high enough and then give an interesting guideline for the study of slip in FE simulations.

Introduction.

We study the phenomenon of wall slip that appears in the simulation of viscoelastic fluids obeying the upper convected Maxwell (UCM) model and the corotational Maxwell (COM) model, in the case of a flow through a 4 : 1 contraction (cf. [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF]). In the case of the UCM model, this phenomenon appears with the use of the new finite element method (FEM) that we have developed in [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF], so-called θ-MSUPG method. We study the influence of the splitting and the modified upwind term in the occurrence of this slip. To do so, we give a comparative study between the θ-MSUPG method and the θ-SUPG method, which uses a classical unmodified SUPG decentering term. To study the special case θ = 1, without the use of splitting, we use a new FE satisfying the inf-sup Tensor-Velocity conditions (P 1 ⊕ b 4 continuous FE approximation of the tensors). The classical Galerkin method and the classical SUPG method are included in this study. For the UCM model, we notice that the slip phenomenon does not appear with the SUPG method, which could explain why we do not observe the slip phenomenon in the literature. On the other hand, for the UCM model, it appears numerically that, for the studied FEMs, the upwinding terms used in the θ-MSUPG method are necessary to obtain a slip phenomenon and that the use of splitting of the θ-MSUPG method can effectively stabilize this phenomenon. We give a mesh convergence study for the θ-MSUPG method in the case of a flow in a 4 :1 contraction for the UCM model. We observe a mesh convergence for the θ-MSUPG method, with convergence of the approximate solution, but, without more theoretical knowledge about the continuous problem (which remains a very difficult problem), one cannot conclude to the presence of a numerical artifact or to the existence of an approximate solution which would converge to a solution of the continuous problem. Note that in [START_REF] Trebotich | Toward a solution to the high Weissenberg number problem PAMM[END_REF], in the case of a non zero Reynolds number (Re = 1), the author observes an unresolved structure in the solution downstream from the reentrant corner that may or may not be a numerical artifact. It appears numerically that, for the FEMs studied, the modified upwind term is necessary to obtain a wall slip phenomenon and that the use of splitting can effectively stabilize this phenomenon (see Fig. 3. 19 & Tab. 5). This could explain why this wall slip is not mentioned in the literature review. In this regard, it is difficult to provide a methodological comparison with the literature. This slip does not appear in the simulations based on the principle of the method introduced in [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF] so called the log-conformation method, such as the 2D simulations performed for the Oldroyd model, with a non-zero Newtonian viscosity (see [START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF], [START_REF] Saramito | On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids[END_REF]), or the Leonov model (see [START_REF] Kwon | Finite element analysis of planar 4 :1 contraction flow with the tensor-logarithmic formulation of differential constitutive equations[END_REF]) or such as the simulation of 3D flow in an abrupt contraction for a fluid governed by the Oldroyd or the Maxwell model (see [START_REF] Alves | Application of the log-conformation tensor approach to three dimensional viscoelastic flows[END_REF]). However, with respect to the high Weissenberg number problem, it should be noted that, according to [START_REF] Alves | Application of the log-conformation tensor approach to three dimensional viscoelastic flows[END_REF], no significant improvements appear, with the log-conformation method, in the case of the upper convected Maxwell problem in 3D simulations (see [START_REF] Alves | Application of the log-conformation tensor approach to three dimensional viscoelastic flows[END_REF]). Furthermore, it appears in our simulations, that the determinant of the conformation tensor seems to tend towards 0 at the downstream wall : one of its positive eigenvalues could become zero (see Tab. 4). The method of [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF] requires that this determinant is not zero on the interior of domain (positif in fact). This could explain why the log-conformation method could not take into account the slip at the boundaries.

We note that the situation studied in this paper is different, for example, from that of Poiseuille flows with wall slip, where a slip law is used. In our situation, the problem admits a solution, presenting a slip (or rather a quasi-slip) at the walls of the channel, characterized by a discontinuity of the velocity at each wall and there is no slip law.

A recent work shows the existence of weak theoretical solutions with wall slip for the UCM fluid flow problem through a channel and brings an element of confirmation that slip is not a simulation artifact (see [START_REF] Sandri | Existence and numerical stability study of weak solutions of a Maxwell fluid in a channel[END_REF]). This leads us to make the hypothesis that this type of weak solution would exist in the case of a 4 :1 contraction, without excluding the possible existence of an associated strong solution, with a corner singularity (cf. the studies in [START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF], [START_REF] Lipscomb | Implications of boundary singularities in complex geometries[END_REF]).

For the COM model and a flow in a channel, a wall slip appears for all FEM used (with or without upwind terms, including the usual Galerkin method) from W e 0.66, where W e is the Weissenberg number. This value of W e 0.66 is close to a critical and predictable value of W e c 0.6438. For the value W e = 0.66, we give a mesh convergence study for the θ-MSUPG method and we also observe a mesh convergence in this case.

We give some details on the θ-MSUPG method and on the comparative study of the article. A particularity of the θ-MSUPG method, is to take into account the purely viscoelastic case, the so-called Maxwell model, important in practice. This method presents stability properties which allows to obtain convergence results for a linearized version of the Maxwell model with a priori error estimates (see [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF], [START_REF] Sandri | Sur une méthode d'éléments finis pour les écoulements de polymères[END_REF]). As far as we know, this is the first result of this kind. This method is characterized by the use of a splitting parameter θ ∈]0, 1] and an upwind term λδB(u, τ ) in the test functions which differs in particular from the λδu.∇τ term also used for solving this type of problem (SUPG method (Streamline Upwind Petrov Galerkin)). For 0 < θ < 1, δ = 1-θ θ , the θ-MSUPG method allows to relax the inf-sup velocity-tensor condition given and studied in [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF] for the study of the convergence of the Finite Element (FE) space introduced in [START_REF] Marchal | A new finite element for calculating viscoelastic flow[END_REF]. In this paper, we study the influence of the splitting parameter θ and the influence of the choice of the upwind term on the occurrence of this slip. For that, we use a FE inspired from [START_REF] Ruas | An optimal three-field finite element approximation of the Stokes system with continuous extra stresses[END_REF], satisfying the inf-sup velocity-tensor condition given in [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF]. This FE allows to solve the approximate problem with splitting and also without splitting (case corresponding to θ = 1). We check numerically that, for the UCM model, the use of the λδB(u, τ ) upwind term is necessary to obtain a slip at the downstream wall and that the use of the splitting can effectively stabilize the flow with wall slip. On another hand this slip does not appear with the use of the λδu.∇τ upwind term. We give a heuristic explanation of the extra-stress tensor values found at the downstream wall of the flow in the presence of wall slip for the COM model (see Eqn. (3.6) & Tab. 2) and for the UCM model (see (3.10) & Tab. 4). We also study the influence of the reentrant corner and the use of splitting on the occurrence of wall slip in the COM model case.

We consider the following non-dimensional Oldroyd problem (cf. [START_REF] Bird | Dynamics of polymeric liquids[END_REF]) on a domain Ω ⊂ IR 2 with polygonal boundary Γ (see Fig. 3.12, the geometry of the contraction 4 :1) :

(P)                    σ + λB(u, σ) = 2αd(u) in Ω, -∇.σ -2(1 -α)∇.d(u) + ∇p = f in Ω, ∇.u = 0 in Ω, u = u 0 on Γ 1 , u.n = 0 on Γ 2 = Γ\Γ 1 , [(σ + 2(1 -α)d(u) -pI).n].t = 0 on Γ 2 , σ = σ 0 on Γ -= {x ∈ Γ, λn.u 0 (x) < 0},
where B is given by the objective derivative 4 where Γ -is the inflow section of the boundary (for λ = 0, we then have Γ -= ∅), n is the outward unit normal vector on Γ, t the unit tangent, and n.u 0 the scalar product of n and u 0 . We denote by Γ 1 a subset, with positive measure, of the boundary Γ. Let ds be the measure on the boundary Γ, we assume that the compatibility conditions Γ 1 u 0 .n ds = 0 and u 0 .n = 0 on Γ 2 are satisfied (it is not necessary that u 0 has zero divergence). We assume that the extra-stress tensors σ and σ 0 are symmetric, u is the velocity vector and p is the pressure. We shall assume that α, λ and a are parameters such that 0 < α ≤ 1, λ ≥ 0 and a ∈ [-1, 1].

B(u, σ) = u.∇σ + σω(u) -ω(u)σ -a(d(u)σ + σd(u)), with f ∈ [L 2 (Ω)] 2 , u 0 ∈ [H 1 (Ω)] 2 , σ 0 ∈ [L 2 (Γ -)]
We have denoted for tensors σ, ∇.σ = σ ij,j the divergence of σ, with the Einstein summation convention, σ : τ the sum

σ ij τ ij , ∇.u = u i,i the divergence of u, u.∇σ = u k σ ij,k , ∇p = p, i the gradient of p, d(u) = 1 2 (u i,j + u j,i
) the rate of strain tensor of u and ω(u) = 1 2 (u i,j -u j,i ) the vorticity tensor of u. Other boundary conditions such as u 2 = 0 or the absence of condition on u at the channel exit are also considered in Section 3 ((BC1) or (BC2) conditions).

In the sequel, we denote by ( , ) the scalar product on [L 2 (Ω)] n , n ∈ IN , and we denote by | . | the corresponding norm. For a section of Γ we denote by , the scalar product σ, τ = σ : τ |u 0 .n|ds.

For α = 1 and a = 0, this problem is the corotational Maxwell model (COM Model) and for α = 1 and a = 1, this problem is the upper convected Maxwell model (UCM model).

If we denote η p the polymeric viscosity and η s the solvent viscosity, α is then the viscoelastic fraction of the viscosity : α = η p η p + η s and (1 -α) the Newtonian fraction of the viscosity. The case α = 1 or α close to 1 is important because in many cases the solvent is, for all practical purposes, absent.

In this work, we present numerical results for the method of splitting, with splitting parameter θ ∈]0, 1], described in [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF] for the FE (finite element) approximation of viscoelastic fluids obeying the upper convected Maxwell model. As EVSS (see [START_REF] Rajagopalan | Finite element methods for calculation of steady viscoelastic flow using constitutive equations with a Newtonian viscosity[END_REF]) or DEVSS method (see [START_REF] Fortin | Numerical analysis of the modified EVSS method[END_REF]) (Discrete -Elastic Viscous Split Stress), this method allows to relax the inf-sup tensor-velocity (ISTV) relating tensor and velocity which appears for the Maxwell problem (see [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF], [START_REF] Marchal | A new finite element for calculating viscoelastic flow[END_REF]). This method allows to obtain convergence results for a linearized version of the Maxwell model with a priori error estimates (see [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF], [START_REF] Sandri | Sur une méthode d'éléments finis pour les écoulements de polymères[END_REF]). As far as we know, this is the first result of this kind.

For the FE approximation we use triangular elements with the Taylor-Hood element for the approximation of the velocity-pressure. For the extra-stress-tensor we use P 1 continuous approximation or [P 1 ⊕ b 4 ] continuous FE approximation (we add 3 bubbles functions, see Section 2). For the linear approximation of the extra-stress, the FE problem can be ill-posed for θ = 1, then, the addition of bubbles helps to remedy this. The [P 1 ⊕ b 4 ] 4 FE space is introduced in [START_REF] Ruas | An optimal three-field finite element approximation of the Stokes system with continuous extra stresses[END_REF], with the use of the Crouzeix-Raviart FE for the approximation of (u, p), with discontinuous approximation of the pressure, in place of the Taylor-Hood FE used in this paper.

The outline of the paper is then as follows : in Section 2 we present the FEM (Finite element method) used in this paper, numerical results are described in Section 3. In Section 3.1, we consider a fluid flowing in a rectangle and in Section 3.2, we describe numerical results obtained for the abrupt 4 to 1 contraction. We give also an heuristic explanation of the values of extrastress tensor found at the downstream wall of the flow.

The FEM approximation of (P)

.

Let T = {τ, τ ∈ [H 1 (Ω)] 4 , τ symmetric}, the tensor space, X = {u ∈ [H 1 (Ω)] 2 , u = 0 on Γ 1 , u.n = 0 on Γ 2 } and M = L 2 0 (Ω) = {q ∈ L 2 (Ω), (1, q) = 0} the space of pressures
with zero mean value or M = L 2 (Ω), used, for example, for Neumann outflow boundary, as in Section 3.1. We remark that for the Stokes Problem obtained with λ = 0 in (P), the natural choice for the tensor space T is [L 2 (Ω)] 4 s (s for symmetric). In fact the choice of the space T above (with more regular functions) is made only in view of the discretization with the use of continuous FE. We denote by {T h × X h × M h ⊂ T × X × M } h>0 a family of FE spaces where h is a discretization parameter (in practice, h will be the size of the mesh). With this choice of T , B is then well defined on X × T . As in the approximation of the Stokes problem we shall assume that the following velocity-pressure condition holds :

(2.1) inf

q∈M h sup v∈X h (∇.v, q) |d(v)||q| ≥ β > 0, with β independent of h.
We recall the velocity-stress inf-sup condition on the kernel V h hereunder, introduced in [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF] for the study of the FE developed in [START_REF] Marchal | A new finite element for calculating viscoelastic flow[END_REF], in view to treat the Maxwell problem :

inf v∈V h sup τ ∈T h (τ, d(v)) |d(v)||τ | ≥ β > 0, with β independent of h, (2.2) V h = {v h ∈ X h , (∇.v h , q h ) = 0, ∀q h ∈ M h }.
Now, we give the variational formulation (Q) of problem (P) considered in this paper. Some details to obtain this formulation are given in [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF] and also in [START_REF] Sandri | Sur une méthode d'éléments finis pour les écoulements de polymères[END_REF] where non-homogeneous boundaries conditions and a δ-dependent mesh version are studied (δ is an upwinding parameter) and also in [START_REF] Sandri | Numerical study of a new finite element method for the approximation of viscoelastic fluid flow problems[END_REF] for a numerical study of this method. The constitutive equation is discretized with a modified SUPG (Streamline Upwind Petrov Galerkin) type method denoted MSUPG (modified SUPG) method. We have chosen an upwinding of the kind τ + δλB(u, τ ), with δ > 0. In comparison with the SUPG method, the modification consists to replace the derivative u.∇τ by B(u, τ ), i.e. by the stationary corotational invariant derivative if a = 0 and by the upper convected invariant derivative if a = 1. Let µ ≥ 0 be a parameter, the problem, obtained by introducing a fraction 1 -θ (θ ∈]0, 1]) of the constitutive equation in the momentum equation, is then (see [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF] and [START_REF] Sandri | Sur une méthode d'éléments finis pour les écoulements de polymères[END_REF]) : Find (σ, u, p) ∈ T × (u 0 + X) × M such that ∀(τ, v, q) ∈ T × X × M we have :

(Q)                (σ + λB(u, σ), τ + δλB(u, τ )) + λ(1 + δ) σ, τ Γ - = 2α(d(u), τ + δλB(u, τ )) + λ(1 + δ) σ 0 , τ Γ -, θ(σ, d(v)) + (1 -θ)(2αd(u) -λB(u, σ), d(v)) +2(1 -α)(d(u), d(v)) -(p, ∇.v) + µ(∇.u, ∇.v) = (f, v), (∇.u, q) = 0,
The Dirichlet boundary conditions for σ on Γ -allows to take into account the hyperbolic character of the constitutive equation (see the use of this condition in [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF]), otherwise the constitutive equation is ill-posed. This condition is satisfied in a weak sense in Problem (Q) (see [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF]). The boundary term of the momentum equation is given by the membership u ∈ (u 0 + X). The stabilization term µ(∇.u, ∇.v) keeps the method consistent, because ∇.u = 0. In this paper we consider the corresponding approximation of Problem (Q) so-called θ-MSUPG method : θ-splitting with MSUPG method

Find (σ h , u h , p h ) ∈ T h × (u 0 + X h ) × M h such that ∀(τ, v, q) ∈ T h × X h × M h we have : (Q h )                (σ h + λB(u h , σ h ), τ + δλB(u h , τ )) + λ(1 + δ) σ h , τ Γ - = 2α(d(u h ), τ + δλB(u h , τ )) + λ(1 + δ) σ 0 , τ Γ -, θ(σ h , d(v)) + (1 -θ)(2αd(u h ) -λB(u h , σ h ), d(v)) +2(1 -α)(d(u h ), d(v)) -(p h , ∇.v) + µ(∇.u h , ∇.v) = (f, v), (∇.u h , q) = 0.
In [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF] or [START_REF] Sandri | Sur une méthode d'éléments finis pour les écoulements de polymères[END_REF], a key to obtain some results of convergence, including the Maxwell linearized problem, is to select µ > 0 and to pick in the formulation

θ = 1 1 + δ .
In this paper, the iterative method used to solve Problem (Q h ) is the following implicit scheme.

Let (σ n h , u n h , p n h ) ∈ T h × (u 0 + X h ) × M h be given, (σ n+1 h , u n+1 h , p n+1 h ) ∈ T h × (u 0 + X h ) × M h is then given by : (FP)                      (σ n+1 h + λB(u n h , σ n+1 h ), τ + δλB(u n h , τ )) + λ(1 + δ) σ n+1 h , τ Γ - = 2α(d(u n+1 h ), τ + δλB(u n h , τ )) + λ(1 + δ) σ 0 , τ Γ -, θ(σ n+1 h , d(v)) + (1 -θ)(2αd(u n+1 h ) -λB(u n h , σ n+1 h ), d(v)) +2(1 -α + c)(d(u n+1 h ), d(v)) -(p n+1 h , ∇.v) + µ(∇.u n+1 h , ∇.v) = (f, v) + 2c(d(u n h ), d(v)) (∇.u n+1 h , q) = 0, ∀(τ, v, q) ∈ T h × X h × M h .
The c term allows us to improve, in certain cases and in a very significant manner, the convergence of the fixed point algorithm above (this c term is considered in [START_REF] Keunings | Computer Modeling for Polymer Processing[END_REF], a numerical study relating c to the convergence of (FP) is carried out in [START_REF] Sandri | Numerical study of a new finite element method for the approximation of viscoelastic fluid flow problems[END_REF]). In this paper we shall use values such as c = 0, c = 1 4 , c = 1. We thus have a linear system of the form A(U n )U n+1 = F (U n ). A criterion for stopping the fixed point is to check the residual A(U n+1 )U n+1 -F (U n+1 ) or the increment U n+1 -U n (this last criterion is mainly used in this article). In our computations we have verified that both criteria give similar results.

Except for the continuation methods used at the end of the Section 3.2 for the UCM model, in all the numerical tests of the following Section 3, we shall select (σ 0 h , u 0 h , p 0 h ) be the solution of the problem (Q h ) with λ = 0 in the equations. We remark that (σ 0 h , u 0 h , p 0 h ) is then an approximation of the solution of a Stokes problem.

In view of studying the error and the convergence order of the iterative method (FP), we shall use the following energy norm :

(σ, u, p) = [|σ| 2 + |d(u)| 2 + |p| 2 ] 1/2 .
In the iterative procedure we observe after some iterations a convergence of the kind :

e n ≤ Aγ n , with e n = (σ n h , u n h , p n h ) -(σ n-1 h , u n-1 h , p n-1 h ),
where γ > 0 is a estimation of the rate of convergence which will depend on h in the simulations. We give the curves of convergence of the difference (for n ≥ 2) :

e n = (σ n h , u n h , p n h ) -(σ n-1 h , u n-1 h , p n-1 h )
in log-scale (with the natural logarithm log e ). The fixed-point iterative process is stopped when the values of e n reach a plateau phase of the order of the precision of the machine (see Fig. 3.4, e n e -25 10 -13 for 32-bit computing). We also used relative norm in Section 3.2 (see Fig. 3. [START_REF] Sandri | On a FEM method for a linearized version of the Oldroyd Problem[END_REF], with similar results.

In this paper, we use the two following finite element spaces denoted (FE1) and (FE2). The domain Ω is equipped with a conforming triangulation {T h } h>0 of triangles K. The pair (u, p) is approximated using the Taylor-Hood FE which satisfies (2.1) and the tensor σ is approximated with the P 1 tensor continuous approximation or with the P 1 ⊕ b 4 continuous approximation.

For the FE (FE1), we select the [P 1 ] 4 -[P 2 ] 2 -P 1 continuous FE approximation of (σ, u, p), with P 1 tensor continuous approximation :

T h = {τ ∈ T ; τ K ∈ [P 1 (K)] 4 , ∀K ∈ T h }, X h = {v ∈ X; v K ∈ [P 2 (K)] 2 , ∀K ∈ T h }, (F E1) M h = {q ∈ M ∩ C 0 (Ω); q K ∈ P 1 (K), ∀K ∈ T h },
where, for m ≥ 0 : P m (K) = {v : K → IR , v(x, y) = 0≤i+j≤m α ij x i y j }, and for the FE (FE2), we select the [P 1 ⊕ b 4 ] 4 -[P 2 ] 2 -P 1 continuous FE approximation, with the so-called P 1 -bubbles tensor continuous approximation :

T h = {τ ∈ T ; τ K ∈ [P 1 (K) ⊕ b 4 (K)] 4 , ∀K ∈ T h }, X h , M h , (F E2) with b 4 (K) = span{λ 2 1 λ 2 λ 3 , λ 1 λ 2 2 λ 3 , λ 1 λ 2 λ 2 3 }, where {λ i } 1≤i≤3 is a set of barycentric coordinates of K.
The inf-sup tensor-velocity condition (2.2) is generally necessary if we choose θ = 1 in (Q h ). This condition is, a priori, not satisfied for (FE1) and we can show, using the techniques presented in [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF], that this condition is satisfied for (FE2). In this work we compare also the results obtained with the upwind method where the τ + δλu.∇τ (δ > 0) advection term in the test functions is used in place of the τ + δλB(u, σ) term. This gives the following approximation of Problem (P), so-called θ-SUPG method : θ-splitting with SUPG method

Find (σ h , u h , p h ) ∈ T h × (u 0 + X h ) × M h such that ∀(τ, v, q) ∈ T h × X h × M h we have : (Q h )                (σ h + λB(u h , σ h ), τ + δλu.∇τ ) + λ(1 + δ) σ h , τ Γ - = 2α(d(u h ), τ + δλu.∇τ ) + λ(1 + δ) σ 0 , τ Γ -, θ(σ h , d(v)) + (1 -θ)(2αd(u h ) -λB(u h , σ h ), d(v)) +2(1 -α)(d(u h ), d(v)) -(p h , ∇.v) + µ(∇.u h , ∇.v) = (f, v), (∇.u h , q) = 0,
with (FP ) the iterative method corresponding to (FP), to solve (Q h ). For θ = 1, we find a standard SUPG method, and moreover, if δ = 0, we find the usual Galerkin method.

Numerical results.

In the following Section 3, we study the appearance of the slip for a Poiseuille flow in a rectangle or for a fluid flowing in a contraction 4 : 1. We recall some results obtained in [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF] with the [P 1 ] 4 -[P 2 ] 2 -P 1 continuous FE approximation and we compare these results with numerical results using the [P 1 ⊕ b 4 ] 4 -[P 2 ] 2 -P 1 continuous FE approximation. This last FE space allows us to study the influence of θ, the parameter of splitting, on the convergence of the method. We study also the influence of the choice of the upwind term on the wall slip. A mesh convergence study for the θ-MSUPG method is given in Tab. 1 and 2 (COM model, flow in a rectangle) and Tabs. 3 and 4 (UCM model, contraction 4 : 1). At the end of Section 3.2, continuation method is used for the UCM model. This method allows in particular to reach higher Weissenberg number and to improve the rate of convergence of the fixed point iteration (FP) or (FP ). Let us first consider the Maxwell problem in the case of a Poiseuille type flow. We consider a fluid flowing in the domain Ω

1 =] -∞, +∞[×] -1, 1[ satisfying : (P 1 ) σ + λB(u, σ) = 2αd(u), -∇.σ + ∇p = 0, ∇.u = 0, in Ω 1 , u(x, -1) = u(x, 1) = 0, ∀x ∈ IR .
We assume that the Maxwell fluid, obtained for α = 1, is laminar and moving in the +x direction, and governed by a constant pressure gradient in the flow direction : p, 1 = -f < 0, σ depends only on y. Then, we have the following existence result (see, for example, [START_REF] Guillop É | Global existence and one-dimensional non linear stability of shearing motions of viscoelastic fluids of Oldroyd type[END_REF], Eqn. 2.8 and Proposition 2.1) :

• Case 3.1. Let α = 1 and a = 1 (upper convected Maxwell model). There exists a unique C 1 Poiseuille flow solution of (P 1 ) with u given by

(3.1) u 1 (x, y) = f 2 (1 -y 2 ), u 2 (x, y) = 0. • Case 3.2. Let α = 1 and a = 0 (corotational Maxwell model). If f ≤ f λ , with f λ = 1 2λ
, then there exists a unique C 1 Poiseuille flow solution of (P 1 ) with u given by

(3.2) u 1 (x, y) = y 1 -2f s 1 + 1 -4f 2 λ 2 s 2 ds, u 2 (x, y) = 0.
For λ given, let us search the value of f in the Poiseuille flow of the Maxwell fluid given in the cases 3.1 and 3.2 such that the volume flow rate is 16 3 i.e. the same as in our simulations, then f has to satisfy (we set f λ = +∞ for a = 1) :

∃f ≤ f λ such that 1 -1 u 1 (x, y)dy = 1 -1 4(1-y 2 )dy = 16 3
, with u 1 given by (3.1) or (3.

2). (3.3)

We denote by f the solution of (3.3) and by u the corresponding Poiseuille flow given by (3.1) or (3.2) with f = f . For the UCM model (a = 1), we have f = 8, and for the COM model (a = 0), a numerical study shows that this problem has a unique solution f if and only if λ satisfies (see [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF]) : In our simulation, the slip appears when λ > ∼ λ c . Then, in this case Eq. (3.2) has no solution satisfying the volume flow rate condition.

In [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF], we have numerically studied the case of a flow through a domain Ω denoted also by Q in this section and given by the rectangle Ω = Q = [0, 1] × [0, 8], with Q corresponding to the downstream section of the contraction 4 to 1 (in view of the symmetry of the flow we use the domain Q in place of the domain [-1, 1] × [0, 8]). We have chosen Newtonian Poiseuille flow at the inflow boundary and Neumann condition at the outflow boundary. The problem considered in Q is the following :

(P 2 ) σ + λB(u, σ) = 2αd(u), -∇.σ -2(1 -α)∇.d(u) + ∇p = 0, ∇.u = 0 in Q
with the following boundary conditions :

• at the inflow boundary x = 0 and y ∈ [0, 1] : u 1 (0, y) = 4(1 -y 2 ), u 2 (0, y) = 0 (Poiseuille flow of a Stokes Fluids) and σ = 2αd(u),

• at the outflow boundary x = 8 and y ∈ [0, 1] : either (σ tot .n).n = 0 and u 2 = 0 (see [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF]) denoted (BC1) boundary conditions (fully developed profile for a viscous fluid, see also [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF], Fig. 4) or σ tot .n = 0 denoted (BC2) boundary conditions (no condition is imposed on the velocity field), where

σ tot = (σ + 2(1 -α)d(u) -pI).
We present results obtained with the (BC2) condition.

A comparison with the (BC1) condition is given in Figs. 3.9-3.10. Except at the exit of the channel, computations show that both conditions give similar solutions. For these conditions the space for the pressure is L 2 (Ω).

• at the wall x ∈ [0, 8] and y = 1 :

u 1 (x, 1) = u 2 (x, 1) = 0,
• at the boundary x ∈ [0, 8] and y = 0 (axis of symmetry) :

u 2 = 0, (σ tot .n).t = 0.
The boundary conditions (BC1) are those used in [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF] in the contraction 4 to 1 case. It is also possible to choose u 1 (8, y) = 4(1 -y 2 ), u 2 (8, y) = 0 (Poiseuille flow of a Stokes Fluids) at the outflow boundary, with similar numerical results. With these inflow boundary conditions, we have the same volume rate of flow than in the contraction case. We notice that, with these inflow boundary conditions, the Problem (P) appears as a perturbation of the plane Poiseuille flow. Other inflow boundary conditions such as σ = 0 give similar numerical results. We remark that, in the UCM case the solution of the Poiseuille flow given by u 1 (0, y) = 4(1 -y 2 ) is not the solution of (P), this is due to the choice of the inflow boundary condition σ = 2αd(u). The Weissenberg Number is defined from the value of ∂u 1 ∂y at the outlet by W e = λ| ∂u 1 ∂y (1, 8)| = 8λ, where u is chosen as u 1 (8, y) = 4(1 -y 2 ), u 2 (8, y) = 0 (this corresponds to a Poiseuille flow of a Stokes fluids as outflow boundary condition).

The reliability of the methods was tested on the (P 2 ) problem with boundary conditions such that the solution is the Poiseuille flow given by (3.1) or (3.2) (see also Figs. 3.2, 3.15, 3.27, the cross-sections of u 1 at the middle of the channel or at the middle of the downstream channel of the contraction). The tests give good results. In particular, we have verified that in the case of a P 2 approximation of the constraints, all the methods considered in this paper give the exact solution (3.1) (in this case the exact solution is in the FE space).

The problem (P) is discretized by problem (Q h ) or problem (Q h ), with the Finite Element spaces (FE1) or (FE2). For the computations, we use the triangulation T h given by the mesh M 1 (992 triangles and 529 nodes) (see Fig. 3.1, first graph) which corresponds to the mesh M 2 (see Fig. 3.13) used for the downstream section of the contraction 4 to 1. We also have used successive refined meshes at the wall y = 1 of the initial mesh

M 1 = M 1 1 (see Fig. 3.1, second graph) : meshes M 2 1 , M 3 1 , M 4
1 until mesh M 5 1 (9668 triangles and 5368 nodes). In Table 1, we have reported for each mesh M i 1 , the number of triangles and the number N i of nodes, as well the size (h i , a i ) of the triangle located at the wall y = 1 with vertices (0, 1), ( h i 2 , 1 -a i ), (h i , 1). We also have used the mesh M 3L 1 , which is a more uniform refinement of M 3 1 (5888 triangles and 3095 nodes), this slightly improves the convergence rate of the fixed point method in its final phase (see the curves and the rates of convergence in Fig. 3.6 and Tab. 1).

We denote by (σ Q , u Q , p Q ) the solution of Problem (P 2 ) and by (σ h , u h , p h ) the approximate solution of Problem (P 2 ) given by (Q h ) or by (Q h ). In [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF], the numerical tests show that, for a = 0, a slip appears at the wall of Q when W e becomes greater than k c . We have made simulations around this value.

Let us consider the case a = 0 and W e = 0.64. For these values the solution of (3.3) is then given by f 6.235. Let u = u, where u is given by (3.2) with f = f = 6.235, we find that

|u h -u | ≤ 0.0075 in [2, 6] × [0, 1] (see Fig. 3.2, first graph).
In this way we see that u is stable for u Q .

1 ). We observe a convergence of the FEM solution with the refinement of the meshes towards a regular solution with a smoothing of irregularities. Then, more generally, for the COM model and λ > ∼ λ c , i .e. W e > ∼ k c it appears a slip of the velocity at the wall. For λ λ c , the gradient of pressure p h,1 has a critical value which is close to -1 2λ c . Then, when λ increases, |p h,1 | falls (see Fig. 3.2, second graph), and it appears a slip of the velocity at the wall. These results are valid with the use of the upwind term λδB(u, τ ), they are also valid with the use the upwind term λδu.∇τ of the SUPG method. In Fig. 3.2, the mean of p h,1 is given on the rectangle [START_REF] Bird | Dynamics of polymeric liquids[END_REF][START_REF] Fortin | Numerical analysis of the modified EVSS method[END_REF] × [0, 0.7], where the flow is close to its corresponding Poiseuille flow and where p h,1 is close to a constant.

We give in Fig. 3.4, the curves of convergence of the fixed point iteration method (FP) for the COM model, obtained with θ-MSUPG method and θ-SUPG method. We use FE (FE1) with (θ, δ) = ( 1011 , 1 10 ) and FE (FE2) with (θ, δ) = ( 10 11 , 1 10 ) or (θ, δ) = (1, 1 10 ) (MSUPG or SUPG method) or (θ, δ) = (1, 0) (Galerkin Method). In this Figure, the finite element space (FE1) (resp. (FE2)) is denoted by P 1 (resp. P 1 ⊕ b 4 ).

For W e = 0.63, the solutions and the rates of convergence are similar (see Fig. 3.4, fisrt graph).

For W e = 0.66, a plateau appears in the curves of convergence. It corresponds to the appearance of a slip at the wall (see Fig. 3.4, second graph). It should be noted that the slip appears even for θ = 1 and δ = 0, so without the use of the upwinding term. Tables 1 and2 give a summary of the mesh convergence study performed on the meshes M i 1 . In this study the θ-MSUPG method, with FE space (FE1), is used to solve the COM model with W e = 0.66. We set L = {(x, y), x ∈ [0, 8], y = 1} the boundary {y = 1} of the domain Ω and L = {(x, y), x ∈ [0, 8], y = 0.99}. Let U i = (σ i , u i , p i ) denote the FEM solution obtained on the mesh M i 1 , i = 1, • • • 5. In Table 1, we have reported for each mesh M i 1 the energy norm of U i and various L 2 norms of U i -U i-1 . The coefficients τ c and τ d give an order of convergence of the FEM, for the consecutive columns of Tab. 1, in function of N = N i and h = h i . For example we observe the following mesh convergence in L 2 (Ω) norm : the L 2 (Ω) norm of the difference of the solution U 2 obtained on the mesh M 2 1 with the solution U 1 obtained on the previous mesh M 1 satisfies U 2 -U 1 L 2 (Ω) = 14.37 (see Tab. 1). Then for two successive meshes we have the estimate :

U i -U i-1 L 2 (Ω) = O(N τ c ) = O(N -1.44 i ) = O(h τ d ) = O(h 1.01 i
). We see in Tab. 1, that the energy norm of U i increases with i, U i = O(N 0.76 i ) : this is due to the fact that the velocity has discontinuities. We also observe mesh convergence on the lines L and L in L 2 norm. θ-MSUPG method, FE space (FE1), (θ, δ, µ, c) = ( 1011 , 1 10 , 2, 0). COM m., W e = 0.66, (BC2). (e) L 2 (Ω) norm between solution on the mesh M i 1 and on the previous mesh M i-1 1 . (f ) Idem for L 2 (L) norm or the L 2 (L ) norm for the velocity (the velocity is nul on L). Now we give a link between the Poiseuille flow and the numerical solution when a slip appears. For the Poiseuille flow solution of (P 1 ), the constitutive law gives for the COM model, using the fact that u is independent of x (i.e. we have u 1 (x, y) = u 1 (0, y), ∀x ∈ IR ) (see [START_REF] Guillop É | Global existence and one-dimensional non linear stability of shearing motions of viscoelastic fluids of Oldroyd type[END_REF], Eqn. (2.6)) :

M esh M i 1 M 1 1 M 2 1 M 3 1 (M 3L 1 ) M 4 1 M 5
(f ) u i 1 -u i-1 1 L 2 (L ) - 1.75 0.5 0.1125 0.0495 -2.42 1.76 (f ) u i 2 -u i-1 2 L 2 (L ) - 0.
(3.5) σ 11 (x, y) = λu 2 1,2 (0, y) 1 + λ 2 u 2 1,2 (0, y) , σ 12 (x, y) = u 1,2 (0, y) 1 + λ 2 u 2 1,2 (0, y) , σ 22 (x, y) = - λu 2 1,2 (0, y) 1 + λ 2 u 2 1,2 (0, y) .
Numerically, the value of u 1,2 (0, 1) tends to -∞ and the velocity u tends to the constant (8/3, 0) ((8/3) is the fluid flow rate) (see Fig. 3.3). If we set u 1,2 (0, 1) = -∞ at the wall, formally, the relations above give for the COM model :

σ 11 (x, 1) 1/λ, σ 12 (x, 1) 0, σ 22 (x, 1) -1/λ. (3.6)
We recover numerically these values and then the constitutive equation is satisfied, even at the wall (see Tab 2). θ-MSUPG method, FE space (FE1), (θ, δ, µ, c) = ( 1011 , 1 10 , 2, 0). COM m., W e = 0.66, (BC2). A cross-section of σ 11 for x = 5 is given in Fig. 3.7, for FE space (FE1) and for each mesh M i 2 , 1 ≤ i ≤ 4 (see also Fig. 3.5 for FE space (FE2)). We see that the numerical solution tends to 0 except at the wall. Let L = (x, y) ∈ Q, x = 5, y ∈ [0, 1] denote the section of Q. On L, an artifact appears at the wall, with negative values of σ 11 , as in the the Gibbs phenomenon for the approximation of discontinuous functions. However the L 1 norm of the negative part of σ 11 , denoted by σ - 11 , tends towards 0 with the size of the mesh (see Fig. 3.7). More precisely we have, for the mesh M 1 (resp. for the meshes M 2 1 , M The numerical solution stabilizes itself to the Poiseuille flow until W e = 1.6. We see in Fig 3 .8, comparing results with "P 1 -bubbles" tensor approximation, that the splitting improves in a significant manner the ratio of convergence. We recover this property in Section 3.2, for a fluid flowing in a contraction. In the case of the usual Galerkin method obtained for (θ, δ) = (1, 0), (this method requires here the use of the finite element (FE2)), the fixed point converges for W e = 1 to a regular solution. For W e = 1.6 and the values of c = 1 and 4, this fixed point diverges. We observe here the role of the upwind schemes in the convergence of the fixed point, in the case of the UCM model. For W e = 2.4, for the two kinds of boundary conditions (BC1) and (BC2), a wall slip appears for the MSUPG method (see Fig. 3.9), the first component of the velocity is given in Fig. 3.9. 1011 , 1 10 , 2, 1). MSUPG Upwinding term : δλB(u n h , τ ). FE spaces (FE2). Mesh M 1 . In the second graph in Fig. 3.9, we give the curve of convergence. We see that the fixed point converges slowly (after 100000 iterations). The convergence is relatively unstable with slow oscillations at the begin of the process. In the case of the contraction studied in the next section, we have obtained, for some values of W e, a wall slip and convergence of the fixed point (see, for example, Fig. 3.19, the curve of convergence with a plateau and also Fig. 3.24-3.26, with a continuation method). For the SUPG method, the numerical solution stabilizes itself to the Poiseuille flow given by (3.1), with instabilities at the inflow boundary (see Fig. 3.10, for (BC1) and (BC2) conditions, W e = 1.4). To explain these instabilities, we recall that, with the choice of the inflow boundary condition σ = 2d(u), the solution of (P) differs from the Poiseuille flow. For W e = 3, the iterative process converges for the MSUPG method while it diverges for the SUPG method.

M esh M i 1 M 1 1 M 2 1 M 3 1 (M 3L 1 ) M 4 1 M 5
σ 11 + λu.∇σ 11 -λ(2u 1,1 σ 11 + u 1,2 σ 12 + u 1,2 σ 21 ) = 2u 1,1 , σ 12 + λu.∇σ 12 -λ(u 2,1 σ 11 + ∇.uσ 12 + u 1,2 σ 22 ) = u 1,2 + u 2,1 , σ 12 = σ 21 σ 22 + λu.∇σ 22 -λ(u 2,1 σ 12 + u 2,1 σ 21 + 2u 2,2 σ 22 ) = 2u 2,2 .
(3.7)

For an incompressible laminar Poiseuille flow, the equations (3.7) become, using the relations u.∇σ = 0, ∇.u = 0, u 1,1 = 0, u 2 = 0 :

σ 11 -2λu 1,2 σ 12 = 0, σ 12 -λu 1,2 σ 22 = u 1,2 , σ 22 = 0. (3.8)
we obtain

σ 11 = 2λu 2 1,2 , σ 12 = u 1,2 , σ 22 = 0. (3.9)
This solution corresponds to the values taken in the median part of the flow by σ h given by the θ-MSUPG or θ-SUPG method for W e sufficiently small (Ex : W e = 1.6, see Fig. 3.8). On another hand, if we assume that u 1,2 tends to -∞ at the wall, then, the terms in (3.8) given by

-2λu 1,2 σ 12 , u 1,2 (1 + λσ 22 )
are finite, if and only if σ satisfies at the wall :

σ 12 = 0, σ 22 = - 1 λ . (3.10)
We find similar numerical values for the θ-MSUPG a little after the beginning of the wall. For example, for the UCM model with W e = 2.4, we get at the second quarter of the wall (see also 3.2. Numerical study of the upper convected Maxwell Problem in an abrupt contraction.

We consider the test problem of a fluid flowing in a domain Ω given by an abrupt 4 to 1 plane contraction (see Fig. 3.12), with a domain of length 16 delimited by the vertices (0, 0), (16, 0), (16, 1) (8, 1), [START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF][START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF], (0, 4). We set Γ 1 = Γ\{(x, y) ∈ Γ, y = 0}, where Γ is the boundary of Ω, Γ 2 = {(x, y) ∈ Γ, y = 0}, the axis of symmetry of the contraction and Γ -= {(x, y) ∈ Γ, x = 0}, the inflow boundary.

Fig. 3.12. Geometry of the 4 to 1 contraction.

We select the body force as f = 0, the Dirichtlet conditions for u 1 are chosen as follows (see [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF], for example, for the choice of boundary conditions) :

• At the inflow boundary Γ -: u 1 (0, y) = 1 -y 2 /16, y ∈ [0, 4],

• At the outflow boundary : we choose Dirichtlet conditions (DBC) u 1 (16, y) = 4 -4y 2 , u 2 (16, y) = 0, for y ∈ [0, 1], or the condition (BC1) : u 2 (16, y) = 0, for y ∈ [0, 1], or the condition (BC2) where no conditions on the velocity are imposed, as seen in Section 3.1.

• u 1 (x, 4) = 0, for x ∈ [0, 8], u 1 (8, y) = 0, for y ∈ [START_REF] Alves | Application of the log-conformation tensor approach to three dimensional viscoelastic flows[END_REF][START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF], u 1 (x, 1) = 0, for x ∈ [START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF][START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF].

For u 2 , we select u 2 = 0 on Γ 1 , therefore the condition u.n = 0 on Γ 2 leads to :

• u 2 = 0 on Γ. Then, at the downstream wall we have for the shear stress γ = |u 1,2 (16, 1)| = 8 and the Weissenberg number is then usually defined for this flow as W e = λ γ = 8λ. In our calculations, we also have taken into account the (BC1) and (BC2) boundary conditions of the section 3.1 (where u 1 = 0, or where no velocity condition is imposed at the outflow (see Fig. 3.18, second graph). These boundary conditions give very similar results to those obtained with the (DBC) conditions with respect to the occurrence of slip at the downstream wall in the simulations.

For the UCM model, the inflow boundary condition upon σ is chosen as the value of σ(0, y), where σ is the solution of the Poiseuille flow on IR × [-4, 4] (see also (3.9)). This gives :

σ 11 (0, y) = 2λu 2 1,2 (0, y), σ 12 (0, y) = u 1,2 (0, y), σ 22 (0, y) = 0. The computations are made on the triangulation {T h } h>0 given by the mesh M 2 given in Fig. 3.13. The number of faces of the triangulation is 2310 and the number of nodes is 1221. Computations are also performed on refined meshes at the wall y = 1 of the initial mesh M 2 = M 1 2 and meshes M 2 2 (2812 triangles and 1501 nodes), M 3 2 (3766 triangles and 2035 nodes), M 4 2 (5620 triangles and 3075 nodes) and M 5 2 (9262 triangles and 5121 nodes). In Table 3, we have reported for each mesh M i 2 , the number of triangles and the number N = N i of nodes, as well as the size (h, a) = (h i , a i ) of the triangle K c located at the wall y = 1 with vertices (8, 1), (8 + h 2 , 1 -a), (8 + h, 1). ).

An enlargement of the mesh M 2 is given in Fig. 3.14. The coordinates of the triangle K c located at the reentrant corner, with an edge on the wall y = 1, are resp. about (8, 1), (8.125, 0.979), (8.25, 1) (see Fig. 3.14). For a triangle K, let h K be the diameter of K, let ρ K be the inradius of the incircle in K and let r K be the ratio r K = h K /ρ K . To obtain the slip while minimizing the number of unknowns we used relatively flat triangles on the wall. For example, for the meshes M i 2 , we have r K c ∼ 24.2 at the wall. We have also used the more regular mesh M 3 given in the figure 3.14. For this mesh, the size of the smallest sides of the triangles is the same as for the mesh M 4 2 , but with more vertices (9194 against 5620), the coordinates of K c are (8, 1), (8.0078, 1), (8.0041, 0.9968) with a ratio r K c 4.8. This mesh allows to improve the ratio of convergence of the fixed point iteration obtained for the mesh M 4 2 (see Fig. 3.19, second graph). On the other hand, the singularities observed in the graphs of Fig. 3.21 disappear. Fig. 3.14. First Graph : mesh M 2 at the reentrant corner. Second Graph : mesh M 3 , enlargement of meshes M 3 and M 4 2 at the reentrant corner.

In the case of the θ-MSUPG method applied to the UCM model, for this geometry, the singularity at the corner causes the slip at the downstream wall more easily than in the case of the rectangle. This fact makes it easier to study this wall slip. In the case of the COM model, for the MSUPG and SUPG methods, the wall slip appears for smaller values of W e than in the case of the rectangle (see the study of the COM model at the end of this section).

As in Section 3.1, we consider θ-MSUPG method and θ-SUPG method. We use FE (FE1) with (θ, δ) = (5 6 , 1 5 ) and FE (FE2) with (θ, δ) = ( 5 6 , 1 5 ), (1, 1 5 ), (1, 0). Now, we study the UCM case. Let W e = 0.8, we see in Fig. 3.15, in the case of the mesh M 2 , that the splitting improves the convergence rate : for θ < 1 the convergence is faster than for θ = 1. On another hand, for θ < 1, the convergence rate is similar for the P 1 tensor approximation and for the P 1 ⊕b 4 tensor approximation. This convergence rate is also similar for the two cases of upwind, with λδB(u, τ ) or λδu.∇τ . We remark that the fixed point converges also without upwind (δ = 0). All FEM methods give similar solutions (see the cross-sections of u 1 in Fig. 3.15, second graph, for different FEM methods), including the usual Galerkin method (θ = 1, δ = 0). This shows that simulations with upwind schemes appears reliable compared to the Galerkin method. We give in Fig. 3.16 the graph of u 1 and σ 11 obtained with θ-SUPG method and FE space (FE2) and a comparison with σ 11 obtained with FE (FE1). For σ 11 , we note that the value of the peak at the reentrant corner is larger for the "P 1 -bubbles" tensor approximation than for the P 1 -tensor approximation, we have respectively in Fig. 3.16 : max Ω σ 11 = 18.84 and max Ω σ 11 = 39.86.

Before detailing the case W e = 1, we study the occurrence of slip as a function of W e. : -For W e = 0.8, for all FEM used, there is no wall slip at the downstream channel for the mesh M 2 . However, a wall slip appears as early as the refined meshes M 3 2 for the θ-MSUPG method (the curve of convergence on the refined meshes M 4 2 is given in Fig. 3.19, second graph). The convergence of the iterative method is lost for the SUPG method from the refined mesh M 4 2 , with the FE space (FE2) and the parameters θ = 1, δ = 1 5 (this is the High Weissenberg number problem). For this method, from the mesh M 1 2 to M 3 2 , the peak of σ 11 increases from max Ω σ 11 = 48 to max Ω σ 11 = 164.

-For W e = 1, the θ-MSUPG method, with θ = 5 6 and the FE space (FE1), gives a velocity with a wall slip at the downstream channel for the UCM case (see Fig. 3.17 and the curves of convergence of the first Graph of Fig. 3

.19).

There is no slip phenomenon for the other methods, including the classic SUPG method (necessarily with the use of the FE (FE2)).

-For W e = 1.4, the θ-MSUPG method, with θ = 5 6 and the FE space (FE2), a slip phenomenon also appears (see Fig. 3.22, second graph). ). SUPG Upwinding term : δλu n h .∇τ . FE space (FE2) for u 1 and σ 11 and also (FE1) for σ 11 , mesh M 2 .

We now describe in more detail the case W e = 1, including a study of mesh convergence for the FE space (FE1). ). MSUPG Upwinding term : δλB(u n h , τ ). FE space (FE1). Mesh M 2 .

The slip appears in the iterative process of resolution by propagation from the reentrant corner of the geometry (see Fig. 3.18, first graph, iteration 110). In the COM model case, this slip appears everywhere along the downstream wall in the convergence process (see [START_REF] Sandri | Numerical study around the corotational Maxwell model for the viscoelastic fluid flows[END_REF]). The (BC1) or (BC2) boundary conditions give similar results, with a quasi-constant profile at the outflow boundary for the (BC2) condition (see Fig. 3.18, second graph, the first component of the velocity solution obtained with the mesh M 3 2 ). For these conditions, the convergence curves are then very close to the curves obtained for condition (DBC). The tensor σ 11 gives, on the wall y = 1, four oscillations in a thin layer of triangles that are adjacent to this wall (see Fig. 3.17, second graph, see also the cross-sections of σ 11 (mesh M 2 ) in Fig. 3.20, see also Fig. 3.21, first graph). We notice that these oscillations correspond to the geometry of the mesh and to the approximation of very irregular functions with discontinuities. We believe that this is also due to the three small irregularities on u 1 that appear at this boundary (see the cross-sections of u 1 (mesh M 2 ) in Fig. 3.20). We notice that this irregularities disappear on the more and thinner regular mesh M 3 (see Fig. 3.20 (mesh M 3 ) and Fig. 3.21, first graph). In the mesh convergence process, we see that these oscillations disappear little by little on the successive refined meshes M i 2 , i = 2, 3, 4 (see Fig. 3.21, second graph). On the other hand, we also repositioned the two triangles with a common vertex of coordinates (2, 1) or (4, 1) or (6, 1) located on the boundary (see Fig. 3.14), in order to obtain a uniform mesh along the downstream boundary. We see in Fig. 3.20 (dashed curve) that the oscillations on σ 11 disappear significantly. For the θ-MSUPG method with FE space (FE1), we see in Fig. 3.19, second graph, that the rate of convergence increases with the mesh refinement. On the mesh M 4 2 , the convergence process ends after a long plateau (see Fig. 3.19, second graph). To improve the convergence on this mesh, we also use the following continuation method, starting from mesh M 3 2 to M 4 2 : we solve the Maxwell problem on the mesh M 3 2 , then, to solve on the mesh M 4 2 , we choose as the starting iteration (σ 0 h , u 0 h , p 0 h ), the interpolation of this solution to the DoF (degrees of freedom). On the mesh M 5 2 , the convergence error oscillates around a plateau, with a solution that remains similar to that shown in Fig. 3.17, the processus of convergence is accelerated with c = 0 and larger value of c > 1/4, c = 0.371 allows to obtain convergence. With the refinement of mesh, the wall slip phenomenon appears for lower Weissenberg numbers : for the value W e = 0.6 on the mesh M 4 2 and for the value W e = 0.8 on the meshes M 3 2 and M 4 2 . Tables 3 and4 give a summary of the study of mesh convergence on meshes M i 2 , for W e = 1. We set L = {(x, y), x ∈ [START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF][START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF], y = 1} the boundary {y = 1} of the domain Ω and L = {(x, y), x ∈ [START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF][START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF], y = 0.99}. Let U i = (σ i , u i , p i ) denote the FEM solution obtained on the mesh M i 2 , i = 1, • • • , 5. The coefficients τ c and τ d are defined as in Section 3.1 (orders of convergence). For the θ-MSUPG method we observe the following mesh convergence in L 2 (Ω) norm : the norm of the difference of the solution U 2 obtained on the mesh M 2 2 with the solution U 1 obtained on the previous mesh M 1 2 is 8.2826. This norm becomes 3.3677 for the meshes M 2 2 -M 3 2 and 2.0102 for the meshes M 3 2 -M 4 2 . Then for two successive meshes we have the estimation :

U i -U i-1 L 2 (Ω) = O(N τ c ) = O(N -1.38 i ) = O(h τ d ) = O(h 0.85 i
) which is close to the estimation reported in Tab. 1 (COM model). The energy norm of U i increases with i, U i = O(N 0.91 i ). We also observe mesh convergence on the lines L and L in L 2 norm. θ-MSUPG method, FE space (FE1), (θ, δ, µ, c) = ( 56 , 1 5 , 2, 1 4 ). UCM m., W e = 1, (DBC). M esh (e) L 2 (Ω) norm between the solution on the mesh M i 2 and on the previous mesh M i-1 2 , (f ) idem for the L 2 (L) norm or the L 2 (L ) norm for the velocity (the velocity is nul on L).

M 1 2 = M 2 M 2 2 M 3 2 M 4 2 M 5
A cross-section of each component of σ and of the determinant of the conformation tensor c = σ + 1 λ Id is given in Fig. 3.21, second graph, for different meshes. In this Figure, we observe the convergence of the tensor at the wall for 9 ≤ x ≤ 15, y = 1 towards a constant. Now, as in Section 3.1, we study the relationship between the constitutive equation and the equation satisfied at the wall by σ, when the flow exhibits a slip at the downstream wall. In this case, we find an approximation of the values (see Eqn. (3.10)) : σ 12 = 0, σ 22 -1 λ satisfied at the wall (see also Eqn. (3.11)). For example, for the UCM model with W e = 1, we have for the mesh M 3 (see Fig. 3.21, second graph, the graph of σ 11 obtained on this mesh) : These values are given in Tab. 4, for the meshes M i 2 .

the solution obtained with the FE space (FE1). In order to study the issue of numerical artifacts, we study the incompressibility of the numerical solution and the residual of the constitutive equation, computed without the B term. In Figure 3.23, we study the values of ∇.u at the nodes and the residual R of the constitutive equation, computed without the B term, with R = A(U n )U n -F (U n ) where A is computed without terms of upwinding. These values are of the same order as in the case of the Stokes's Problem (for ∇.u) or for simulations with the SUPG Upwinding term. When there is a wall slip, the highest values are at the downstream wall, instead only of the corner singularity. In Figure 3.23, we give the maximum, minimum and average values of |∇.u| and |R|. For the L 1 norm of |∇.u|, the respective values of this norm for the Stokes Problem, the θ-SUPG and the θ-MSUPG methods are respectively Ω |∇.u| = 0.2007, 0.2289, 0.4869. On Ω the domain Ω without the triangles having at least one vertex located on the downstream wall, we have Ω |∇.u| = 0.0185, 0.02, 0.0649. In the sequel we study the convergence of the fixed point iteration (FP) or (FP ), with the following continuation method : we choose an initial W e and we solve the Maxwell problem for this number W e. The obtained solution is then the starting iterate in (FP) or (FP ), with W e = W e + 0.1. For example, in the results of the left graph of Fig. 3.24, we first solve the problem for W e = 0.8 and we plot the convergence curve for this value of W e. Then, choosing as starting iteration of the fixed point, this solution obtained for W e = 0.8, we solve the problem for W e = 0.9, and we plot the convergence curve obtained. Then, according to the same principle, we solve for W e = 1, 1.1, etc. until W e = 2.3, value for which the scheme diverges (plateau with oscillations). This method allows in particular to reach higher Weissenberg numbers and to improve the rate of convergence of the fixed point iteration. For example, in the case of the θ-MSUPG method, with FE space (FE2), we reach W e = 2.3. We give the curves of convergence of the continuation method for the FE methods studied in this section starting from W e = 0.8 for the θ-MSUPG method and from W e = 1 for the θ-SUPG method. For the θ-MSUPG method with θ < 1 and FE space (FE1), the (FP) iteration converges until W e = 2.2 and gives a solution with wall slip for the velocity from W e = 1 (see Fig. 3.24, first graph). Then, for W e = 2.3, the convergence error oscillates around a plateau, with a solution that remains similar to that shown in Fig. 3.17 and the process is stopped. We observe a plateau in the curve of convergence when the slip appears for W e = 1. The process is similar for the FE space (FE2) (see Fig. We now give the corresponding convergence curves for the θ-SUPG method, with an initial Weissenberg number of W e = 1. For θ < 1 and the FE espace (FE1) (resp. (FE2)), we obtain a convergence of the (FP ) iteration from W e = 1 until W e = 1.5 (resp.W e = 1.7), without any slip phenomenon (see Fig. 3.25). Then, the iterative process diverges for W e = 1.6 (resp. W e = 1.8). In both cases, the solutions diverge. We study the case where there is no splitting. For the θ-MSUPG method with θ = 1 and FE space (FE2), we see that an unstable slip appears but only for W e = 1.1 (the graph of the solution is close to that of the solution given in Fig. 3.22, first graph), then the iterative method becomes unstable (see Fig. 3.26, first graph) and the approximate solution diverges. There are no slip phenomenon for the SUPG method (see Fig. 3.26, second graph). For these two cases, the solutions diverge. In conclusion, we see, that, even if the ISTV condition is satisfied, the splitting always remains necessary to obtain a stable wall slip. We give a summary of the above results on the mesh M 2 , extrapolating the results to the other W e numbers for which the calculation has not been performed.

θ-MSUPG method θ < 1, FE space (FE1) : no slip for W e ∈ [0, 0.9], slip for W e ∈ [1, 2.2], plateau for W e = 2. Then, for W e ≥ 0.62 a slip appears at the downstream channel for all the methods, and thus before the critical value W e c 0.6438 encountered in the case of a flow in a rectangle (see Eqn. (3.4)). Therefore, in the contraction case, we see that the critical Weissenberg number for the occurrence of the slip decreases, this seems mainly due to the influence of the reentrant corner.

3. 1 .

 1 Slip for the Maxwell Model (corotational and upper convected model) for a Poiseuille flow.

(3. 4 )

 4 λ ≤ λ c where λ c 0.080475 i .e. W e ≤ k c where k c 0.6438, and for λ = λ c we have f = f λ c =

Fig. 3 .

 3 Fig. 3.1. Mesh M 1 (the size of the mesh of the wall y = 1 is 1 4 ) and a portion of the meshes M 2 1 , M 3 1 , M 4 1 , M 3L 1 (the respective sizes of the meshes at the wall are h = (1 16 , 1 32 , 1 32 , 1 64 , 1 128 ).

Fig. 3 .

 3 Fig. 3.2. COM model, (BC2). First Graph, cross-section of u 1 on the line x = 4, W e = 0.64. Legend :section of u 1h the first component of the approximate velocity. o : Exact Poiseuille flow given by (3.2) with f = 6.235. Second Graph, Gradient -p h,1 of pressure vs. W e. Legend : -× --mean value of -p h,1 on [2, 6] × [0, 0.7], -Critical value of case (3.2) : f c = 1 2λ = 4 W e . -• --Value of f solution of (3.2), (3.3), for W e given. Parameters of the FEM and of the fixed point : (θ, δ, µ, c) = ( 10 11 , 1 10 , 2, 0), MSUPG Upwinding term : δλB(u n h , τ ), FE spaces (FE1). Then, more generally, for the COM model and λ < ∼ λ c , i .e. W e < ∼ k c , the numerical solution stabilizes itself to the solution of the corresponding Poiseuille flow, furthermore p h the approximate pressure satisfies p h,1 -f . Now, let W e = 0.66, for this case Eqn. (3.3) has no solution u given by Eqn. (3.2), we see that a slip appears in the domain Q (see Fig. 3.3, the graph of u 1 obtained with the meshes M 1 , M41 ). We observe a convergence of the FEM solution with the refinement of the meshes towards a regular solution with a smoothing of irregularities.

Fig. 3 .

 3 Fig. 3.3. Cross-section of u 1 on the line x = 4 (approximate velocity) and graph of u 1 (reverse angle). COM model with W e = 0.66, (BC2). Parameters : (θ, δ, µ, c) = ( 10 11 , 1 10 , 2, 0). MSUPG Upwinding term : δλB(u n h , τ ). FE spaces (FE1). Results on mesh M 1 and mesh M 4 1 .

Fig. 3 . 4 .

 34 Fig. 3.4. Curves of convergence. COM model with W e = 0.63 and W e = 0.66, (BC2). Parameters : µ = 2, c = 0, (θ, δ) = ( 10 11 , 1 10 ) and for the FE space (FE2) only : (θ, δ) = (1, 1 10 ) or (θ, δ) = (1, 0). SUPG and MSUPG Upwinding terms. Mesh M 1 (32-bit computing).

  (a) base h × altitude a of the triangle K c located along the wall y = 1. (b) Estimation of the rate of convergence (slope of the curves in Fig.3.6).(c) Approximation of the values of the rate of convergence γ from N (number of nodes).(d) Energy norm or L 2 norm of the FEM solution on mesh M i 1 .

Tab. 2 .

 2 Mean values of the solution at the wall y = 1 or y = 0.99. As in Tab. 1 (see the rows denoted by (f ) ), the terms τ c and τ d are the order of convergence in L 2 ( ) norm, but with = {(x, y), x ∈ [1, 7], y = 1} for the tensor and the pressure and = {(x, y), x ∈ [1, 7], y = 0.99} for the velocity.

Fig. 3 .

 3 Fig. 3.7. Cross-section of σ 11 at x = 5 for meshes M 1 , M 2 1 , M 3 1 , M 4 1 . COM model, W e = 0.66, (BC2). Parameters : (θ, δ, µ, c) = ( 10 11 , 1 10 , 2, 0). Upwinding term : δλB(u n h , τ ). FE spaces (FE1).

  Now, if we consider the UCM model (a = 1) then (3.3) is always satisfied with f = 8 in (3.1).

Fig. 3 . 8 .

 38 Fig. 3.8. Curves of convergence. UCM model, W e = 1.6, (BC2). Parameters : (µ, c) = (2, 1), (θ, δ) = ( 10 11 , 1 10 ) and for the FE space (FE2) only, (θ, δ) = (1, 1 10 ) or (θ, δ) = (1, 0). SUPG and MSUPG Upwinding terms. Mesh M 1 .

Fig. 3 . 9 .

 39 Fig. 3.9. First graph : first component of the velocity (reverse angle). Second graph : Curves of convergence. UCM model with W e = 2.4, (BC1) and (BC2). Parameters : (θ, δ, µ, c) = (1011 , 1 10 , 2, 1). MSUPG Upwinding term : δλB(u n h , τ ). FE spaces (FE2). Mesh M 1 . In the second graph in Fig.3.9, we give the curve of convergence. We see that the fixed point converges slowly (after 100000 iterations). The convergence is relatively unstable with slow oscillations at the begin of the process. In the case of the contraction studied in the next section, we have obtained, for some values of W e, a wall slip and convergence of the fixed point (see, for example, Fig.3.19, the curve of convergence with a plateau and also Fig.3.24-3.26, with a continuation method). For the SUPG method, the numerical solution stabilizes itself to the Poiseuille flow given by (3.1), with instabilities at the inflow boundary (see Fig.3.10, for (BC1) and (BC2) conditions, W e = 1.4). To explain these instabilities, we recall that, with the choice of the inflow boundary condition σ = 2d(u), the solution of (P) differs from the Poiseuille flow. For W e = 3, the iterative process converges for the MSUPG method while it diverges for the SUPG method.

Fig. 3 .

 3 Fig. 3.10. First graph : first component of the velocity (reverse angle). Second graph : component σ 22 of the tensor. UCM model with W e = 1.4, (BC1) and (BC2). Parameters : (θ, δ, µ, c) = ( 10 11 , 1 10 , 2, 1). SUPG Upwinding term : δλu n h ∇τ . FE spaces (FE2). Mesh M 1 . Now we give for the UCM model, as in the COM model, a links with the Poiseuille flow and the numerical solution when a slip appears. The constitutive equations of the Maxwell model

Fig. 3

 3 Fig. 3.11. Graph of the components of σ and det(σ + 1 λ Id) (reverse angle). UCM model with W e = 2.4, (BC2). Parameters : (θ, δ, µ, c) = ( 10 11 , 1 10 , 2, 1). MSUPG Upwinding term : δλB(u n h , τ ). FE spaces (FE2). Mesh M 1 .

Fig. 3 .

 3 Fig. 3.13. Initial mesh M 2 (the size of the mesh of the wall y = 1 is 1 4) and mesh M 4 2 (the size of the mesh along the wall y = 1 is 1 32 ).

Fig. 3 .

 3 Fig. 3.15. UCM, W e = 0.8, (DBC). First graph : curves of convergence. Parameters : (µ, c) = (2, 1 4 ), (θ, δ) = ( 5 6 , 1 5 ) and for the FE space (FE2) only : (θ, δ) = (1, 1 5 ) or (θ, δ) = (1, 0). Second graph : section of u 1 at x = 12 (middle of the downstream channel) for different FEM, the approximate solutions u 1h are close to the Poiseuille flow u 1 (x, y) = 4(1 -y 2 ), mesh M 2 .

Fig. 3 .

 3 Fig. 3.16. Graph of the first component of the velocity (reverse angle) and comparative graphs of the solution σ 11 obtained with the space (F E1) space or the space (F E2). UCM model, W e = 0.8, (DBC). Parameters : (θ, δ, µ, c) = (56 , 1 5 , 2, 1 4 ). SUPG Upwinding term : δλu n h .∇τ . FE space (FE2) for u 1 and σ 11 and also (FE1) for σ 11 , mesh M 2 .

Fig. 3 .

 3 Fig. 3.17. Graphs of u 1 and σ (reverse angle). UCM model, W e = 1 (DBC). Parameters : (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4). MSUPG Upwinding term : δλB(u n h , τ ). FE space (FE1). Mesh M 2 .

Fig. 3 .

 3 Fig. 3.18. First component of the velocity (reverse angle). UCM model, W e = 1. Parameters : (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4 ). MSUPG Upwinding term : δλB(u n h , τ ). FE space (FE1). First graph : iteration number 110 (DBC), mesh M 2 . Second graph : (BC2) boundary conditions, mesh M 3 2 . The appearance of the slip is characterized by a plateau in the curve of convergence (see the fisrt graph of Fig. 3.19). Comparing the rate of convergence obtained for W e = 0.8, we see that the rate of convergence increases slightly with W e.

Fig. 3 . 19 .

 319 Fig. 3.19. Curves of convergence. UCM model (DBC). First graph : W e = 1, mesh M 2 . Parameters : θ = 5 6 or θ = 1 (for the FE space (FE2) only), (δ, µ, c) = ( 1 5 , 2, 1 4 ). Second graph : curves of convergence in relative norm for W e = 1 on meshes M 2 , M 2 2 , M 3 2 , M 4 2 , M 3 and for W e = 0.8 and W e = 0.6 on mesh M 4 2 . θ-MSUPG method with FE space (FE1). Parameters : (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4 ).

Fig. 3 .Fig. 3 .

 33 Fig. 3.20. Cross-sections of u 1 and σ 11 . UCM model, W e = 1 (DBC). Parameters : (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4). MSUPG Upwinding term : δλB(u n h , τ ). FE space (FE1). Mesh M 2 and mesh M 3 .

  ) denotes the FEM solution on the meshM i 2 , L = {(x, y), x ∈ [8, 16], y = 1}, L = {(x, y), x ∈ [[START_REF] Hinch | The flow of an Oldroyd fluid around a sharp corner[END_REF][START_REF] Moreno | Logarithmic conformation reformulation in viscoelastic flow problems approximated by a VMS-type stabilized finite element formulation[END_REF], y = 0.99}.(a) base h × altitude a of the triangle K c located at the re-entrant corner. (b) Estimated rate of convergence γ (slope of the curves in the second graph of the Fig.3.19). (c) Approximation of the value of the rate of convergence γ from N (number of nodes).(d) Energy norm or L 2 norm of the FEM solution on mesh M i 2 .

  (x, 1)dx = -0.2804 0.

Fig. 3 .

 3 Fig. 3.22. First component of the velocity (reverse angle). UCM model (DBC). Parameters : (δ, µ, c) = ( 1 5 , 2, 1 4 ). MSUPG Upwinding term : δλB(u n h , τ ). FE space (FE2). Mesh M 2 . First graph : W e = 1.1, θ = 1. Second graph : W e = 1.4, θ = 5 6 .

Fig. 3 .

 3 Fig. 3.23. First Figure. Isovalues of |∇.u|, for the Stokes Problem (First Graph), for the UCM model with W e = 1.4, with solutions obtained with SUPG Upwinding term δλu n h .∇τ (Second Graph) or MSUPG Upwinding term δλB(u n h , τ ) (Third Graph). Parameters (δ, µ, c) = ( 1 5 , 2, 1 4 ). Second Figure. Isovalues of the residual : R = A(U n )U n -F (U n ) with A computed without upwinding, i.e. δ = 0 in (FP) or (FP ). We take on each triangle K, R = max i |R(i)|, with i ∈ DoF of K for the constitutive equation. (DBC), FE spaces (FE2), mesh M 2 .

  3.24), with convergence and wall slip phenomenon from W e = 1.4 until W e = 2.5. For W e = 2.6, the curve of convergence oscillates around a plateau, with a solution close to that in Fig.3.22, second graph.

Fig. 3 .

 3 Fig. 3.24. Curves of convergence, continuation method. UCM model and W e ≥ 0.8 (DBC). Step size for W e : 0.1. Approximation with θ-MSUPG method. Parameters (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4 ). First graph : FE space (FE1), no slip for W e ∈ [0.8, 0.9], slip for W e ∈ [1, 2.2], plateau for W e = 2.3. Second graph : FE space (FE2), no slip for W e ∈ [0.8, 1.3], slip for W e ∈ [1.4, 2.5], plateau for W e = 2.6, mesh M 2 .

Fig. 3 .

 3 Fig. 3.25. Curves of convergence, continuation method. UCM model, W e ≥ 1 (DBD), step size for W e : 0.1. Approximation with θ-SUPG method. Parameters (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4 ). First Graph : FE space (FE1), convergence for 1 ≤ W e ≤ 1.5. Divergence for W e = 1.6. Second graph : FE space (FE2). Convergence for 1 ≤ W e ≤ 1.7. Divergence for W e = 1.8, mesh M 2 .

Fig. 3 .

 3 Fig. 3.26. Curves of convergence, continuation method. UCM model. W e ≥ 0.8, step size for W e : 0.1. Parameters (θ, δ, µ, c) = ( 5 6 , 1 5 , 2, 1 4 ). FE spaces : [P 1 ⊕ b 4 ] 4 (tensor)-[P 2 ] 2 (velocity)-P 1 (pressure). First Graph : upwinding term : δλB(u n h , τ ) (MSUPG), no slip for W e ∈ [0.8, 1], unstable slip for W e = 1.1, divergence for W e = 1.2. Second Graph : upwinding term : δλu n h ∇τ (SUPG), convergence for W e ∈ [1, 1.7] (no slip), divergence for W e = 1.8, mesh M 2 .

3 .θ = 1 ,

 31 FE space (FE2) : no slip for W e ∈ [0, 1], unstable slip for W e = 1.1, div. for W e = 1.2. θ < 1, FE space (FE2) : no slip for W e ∈ [0, 1.3], slip for W e ∈ [1.4, 2.5], plateau for W e = 2.6. SUPG method θ < 1, FE space (FE1) : for W e ∈ [0, 1.5] : convergence and no slip, for W e = 1.6 : divergence. θ = 1, FE space (FE2) : for W e ∈ [0, 1.7] : convergence and no slip, for W e = 1.8 : divergence. θ < 1, FE space (FE2) : for W e ∈ [0, 1.7] : convergence and no slip, for W e = 1.8 : divergence. Tab. 5. Comparative study of θ-MSUPG and θ-SUPG methods.We now study the influence of the reentrant corner for the COM (corotational Maxwell) model. For the COM model, the inflow boundary condition upon σ is chosen as the value of σ(0, y) given in Eqn.(3.5). For W e = 0.40, all the solutions are similar and there is no slip of the velocity at the downstream wall. For W e = 0.58, we have two kinds of FEM solutions : for θ < 1, at the downstream channel, the solutions are close to the Poiseuille flow of the COM fluid given by Eqn. (3.2) (see Fig.3.27, the cross-sections of u 1 at x = 12, for W e = 0.58). On another hand, for θ = 1, a slip appears at the downstream wall (see Fig.

  3.27), including the usual Galerkin method (θ = 1, δ = 0).

Fig. 3 .

 3 Fig. 3.27. COM model, W e = 0.58 (DBD) for different FEM. First graph : curves of convergence. Parameters : (δ, µ, c) = ( 1 5 , 2, 0), θ = 5 6 and only for the FE space (FE2) : θ = 1. Second graph : section of u 1 at x = 12, for different FEM. We have two kinds of solution u 1h . Mesh M 2 .

  Tab. 3. Study of mesh convergence. U i = (σ i , u i , p i

								2	τ c	τ d
	Number of triangles	2310	2812	3766	5620	9262
	Number of nodes N i	1221	1501	2035	3075	5121
	Number of unknowns	14386	17630	23810	35838	59490
	(a) size h × a			1 4 × 1 48	1 8 × 1 96	1 16 × 1 192	1 32 × 1 384	1 64 × 1 768
	(b) γ			0.9682	0.9932	0.9985	0.9996	0.9999
	(c) γ ∼1-1.15e11×N -4.14 i	0.9807	0.9918	0.9977	0.9996	0.9999
	(d) Energy norm of U i	65.7551	90.5165	126.1757	177.6326	251.1289	0.91 -0.48
	(d) L 2 (Ω) norm of U i 33.7177	29.1622	27.9855	27.4314	27.1067
	(e) U i -U i-1		L 2 (Ω)	-	8.2826	3.3677	2.0102	1.3884	-1.38 0.85
	(f ) σ i 11 -σ i-1 11	L 2 (L)	-	43.1056	15.6655	9.8178	5.8908	-1.52 0.93
	(f ) σ i 12 -σ i-1 12	L 2 (L)	-	1.9241	1.3334	0.6728	0.4511	-1.22 0.73
	(f ) σ i 22 -σ i-1 22	L 2 (L)	-	1.0879	1.8865	0.9477	0.7913	-0.44 0.24
	(f ) u i 1 -u i-1 1		L 2 (L )	-	1.6717	0.2433	0.0547	0.0263	-3.29 2.01
	(f ) u i 2 -u i-1 2		L 2 (L )	-	0.0465	0.0503	0.0284	0.0099	-1.33 0.75
	(f ) p i -p i-1	L 2 (L)	-	6.1230	3.4053	2.4665	1.9437	-0.89 0.54

) and for the FE space (FE2) only : (θ, δ) = (1, 1 5 ) or (θ, δ) = (1, 0). Second graph : section of u 1 at x = 12 (middle of the downstream channel) for different FEM, the approximate solutions u 1h are close to the Poiseuille flow u 1 (x, y) = 4(1 -y 2 ), mesh M 2 .

Tab. 4. Mean values of the solution at the wall. As in Tab 3 (see the rows denoted by (f ) ), the terms τ c and τ d are the order of convergence in L 2 ( ) norm, but with = {(x, y), x ∈ [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF][START_REF] Marchal | A new finite element for calculating viscoelastic flow[END_REF], y = 1} for the tensor and the pressure and = {(x, y), x ∈ [START_REF] Johnson | Finite element methods for linear hyperbolic problems[END_REF][START_REF] Marchal | A new finite element for calculating viscoelastic flow[END_REF], y = 0.99} for the velocity.

The determinant of the conformation tensor seems to tend towards 0 at the downstream wall : one of its positive eigenvalues could become zero (this could correspond to a type change of the equations, but only at the wall).

Let us consider the θ-MSUPG method with the FE space approximation (FE2). For θ = 1, a slip phenomenon with instabilities appears also at the downstream wall for W e = 1.1, but only until the middle part of the downstream wall (see Fig. 3.22, first graph). Then, from W e = 1.2, the error of convergence oscillates around a plateau (see Fig 3.26, first graph, the continuation method case), an unstable slip appears at the beginning of the iterative process (as in Fig. 3.22, first graph) and then the approximate solution diverges. This shows that for the FEM used in this paper, the splitting is not necessary to obtain the beginning of a slip phenomenon but it is still necessary to obtain a slip all along the downstream wall of the flow for various W e number.

For θ < 1, a slip phenomenon also appears at the downstream wall for W e = 1.4 (see Fig. 3.22, second graph). The converged solution u 1 presents an instability at the beginning of the downstream wall. This instability disappears, if we choose as starting iteration (σ 0 h , u 0 h , p 0 h ) be

We give the graph of the first component of the velocity and the graph of σ 11 in Fig. 3.28.

Conclusion.

For the COM model and for all FEM used, including the usual Galerkin method, a phenomenon of wall slip appears from a critical Weissenberg W e c . For the flow in a channel, the apparition of this phenomenon is predictable : it corresponds to a loss of solution in the case of the Poiseuille flow (see Eqn. (3.2), (3.3)). At the wall, the comparisons of the numerical values given by Tab 2. and values given by (3.6) are found to be in good agreement. We observe a mesh convergence for the θ-MSUPG method, with a factor of convergence of the iterative method satisfying γ = 1 -O(N -2.1 ), with N the number of nodes of the mesh. The convergence order of the FEM is estimated by O(N -1.44 ) in L 2 (Ω) norm. For the COM model, we found similar results of phenomenon of wall slip in the case of the 4 :1 contraction. This wall slip phenomenon also appears for the UCM model, but only for the θ-MSUPG method. In the case of the flow in a channel, this wall slip appears for a high Weissenberg number, with a relatively unstable solution. For flow through the 4 :1 contraction with W e = 1 and the FE (FE1), we observe a mesh convergence for the θ-MSUPG method, with slip at the downstream wall. The convergence order of the FEM is estimated by O(N -1.38 ) in L 2 (Ω) norm. If we sufficiently refine the meshes, the slip phenomenon also appears for smaller values of W e = 0.8, 0.6 as well as for W e = 1 with the FE (FE2). It appears that the MSUPG-upstream term and the splitting are necessary to obtain solution with a wall slip. Apart from the fact that the θ-MSUPG method is the only method among those used which allows to obtain convergence results in the linearized version of the Maxwell problem, without more theoretical results : for example, whether the UCM problem admits solutions or not, and without a priori convergence results of the FEM, we cannot conclude whether we have a numerical artifact or not. It is as if the problem admits a strong solution which is repulsive for the fixed point and also admits a solution with slip which is attractive for the fixed point, provided that the FE method used allows the approximation of solutions with slip. However, a recent work shows the existence of weak theoretical solutions with slip for the UCM fluid flow problem through a channel, in addition to the classical strong solutions, and provides an element of confirmation that the wall slip observed in the calculations of this paper is not a simulation artifact (see [START_REF] Sandri | Existence and numerical stability study of weak solutions of a Maxwell fluid in a channel[END_REF]).

Computations. All the computations of this paper are made with the FEM code developed in [START_REF] Sandri | Code d'éléments finis pour la résolution de problèmes d'équations aux dérivées partielles. Première partie, mode d'emploi, exemples et éléments mathématiques[END_REF], [START_REF] Sandri | Code d'éléments finis pour la résolution de problèmes d'équations aux dérivées partielles. Seconde partie, algorithmes et description technique du code[END_REF]. This code allows to take into account FE spaces with bubbles functions and more general tests functions like those used in this paper.