Pierre Saramito
email: pierre.saramito@imag.fr

On a modified non-singular log-conformation formulation for Johnson-Segalman viscoelastic fluids

Keywords: Johnson-Segalman viscoelastic fluid, matrix logarithm, Newton method, incompressible finite elements, adaptive mesh, lid-driven cavity

A modified log-conformation formulation of viscoelastic fluid flows is presented in this paper. This new formulation is non-singular for vanishing Weissenberg numbers and allows a direct steady numerical resolution by a Newton method. Moreover, an exact computation of all the terms of the linearized problem is provided. The use of an exact divergence-free finite element method for velocity-pressure approximation and a discontinuous Galerkin upwinding treatment for stresses leads to a robust discretization. A demonstration is provided by the computation of steady solutions at high Weissenberg numbers for the difficult benchmark case of the lid driven cavity flow. Numerical results are in good agreement, qualitatively with experiment measurements on real viscoelastic flows, and quantitatively with computations performed by others authors. The numerical algorithm is both robust and very efficient, as it requires a low mesh-invariant number of linear systems resolution to obtain solutions at high Weissenberg number. An adaptive mesh procedure is also presented: it alows representing accurately both boundary layers and the main and secondary vorties.

Introduction

The Johnson-Segalman [START_REF] Johnson | A model for viscoelastic fluid behavior which allows nonaffine deformation[END_REF] model is considered here, i.e. the upper convective derivative in the Oldroyd-B [START_REF] Oldroyd | On the formulation of rheological equations of states[END_REF] model is replaced by a mixed Gordon-Schowalter [START_REF] Gordon | Anisotropic fluid theory: a different approach to the dumbbell theory of dilute polymer solutions[END_REF] derivative with a parameter a ∈ [-1, 1]:

D a τ Dt = ∂τ ∂t + (u.∇)τ -τ g a (u) T -g a (u) τ , (1)
for all symmetric tensor τ , where u is the velocity field and

g a (u) = 1 + a 2 ∇u - 1 -a 2 ∇u T (2)
denotes a generalized gradient, with the convention ∇u = ∂u i ∂x j 1 i,j d and d = 2, 3. Here, a ∈ [-1, 1] denotes the parameter of the tensor derivative. Remark that g a (u) = W (u) + aD(u) where D(u) = (∇u + ∇u T )/2 and W (u) = (∇u -∇u T )/2 are the symmetric and skew-symmetric parts of the velocity gradient, respectively. When a = 1, then g 1 (u) = ∇u and the Gordon-Schowalter derivative coincides with the usual upper-convected tensor derivative. When a = 0, we obtain the corotational derivative and when a = -1, the lower convected derivative. Problems involving such tensor derivatives appear in non-Newtonian viscoelastic polymer melt flow problems (see e.g. [START_REF] Saramito | Numerical simulation of viscoelastic fluid flows using incompressible finite element method and a θ-method[END_REF][START_REF] Saramito | Méthodes numériques en fluides complexes : théorie et algorithmes[END_REF]), in turbulence modeling with the R ij -ǫ Reynolds stress turbulence models, for liquid crystals [START_REF] Fielding | Nonlinear dynamics and rheology of active fluids: simulations in two dimensions[END_REF], fiber suspension [START_REF] Lipscomb | The flow of fiber suspensions in complex geometries[END_REF] or active fluids [START_REF] Fielding | Nonlinear dynamics and rheology of active fluids: simulations in two dimensions[END_REF], where a is related to the aspect ratio of the particles in suspension. All these models introduce a relaxation time that multiplies the tensor derivative. The corresponding dimensionless quantity is called the Weissenberg number.

The numerical computation of viscoelastic flows in non-trivial geometries has been founded as a very challenging enterprise. The failure of numerical methods when the Weissenberg number becomes large is known as the high Weissenberg number problem. In 1986, Keunings [START_REF] Keunings | On the high Weissenberg number problem[END_REF] observed that the maximal Weissenberg number reached by all algorithms was mesh-dependent: he deduced that this failure was due to a wrong numerical methodology. In 1987, Marchal and Crochet [START_REF] Marchal | A new mixed finite element for calculating viscoelastic flow[END_REF] presented numerical simulations at high Weissenberg numbers obtained with an improved numerical approach and pointed out the need of a mixed finite element methods and of an upwind treatment of the stress transport terms. This was the starting point of many numerical computations of viscoelastic flows at higher Weissenberg numbers (see e.g. [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF][START_REF] Saramito | Efficient simulation of nonlinear viscoelastic fluid flows[END_REF] and [START_REF] Owens | Computational rheology. Imperial college press[END_REF] for a review of this period). In 2000, in a review paper, Keunings [START_REF] Keunings | On the high Weissenberg number problem[END_REF] observed the progress accomplished, but pointed out that the maximal Weissenberg number reached in most reported simulations was still clearly decreasing with mesh refinement. He also concluded on the crucial role of benchmark flow problems for future works.

In 2004, Fattal and Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] observed that some numerical instabilities are caused by the failure of polynomial functions to approximate accurately the exponential growth of the stress tensor, due to the presence of the deformation as a source term in the tensor transport equation. The deformation source term takes its origin in the two last terms on the right-hand side of the tensor derivative [START_REF] Afonso | The finite-volume method in computational rheology[END_REF]. The solution proposed by these authors was a change of unknown that scale logarithmically with the stress tensor. Exploiting the fact that the conformation tensor is symmetric positive definite, the stress transport equation was rewritten as an equation for the matrix logarithm of the conformation tensor (the so-called log-conformation formulation).

Numerical experiments for the driven cavity benchmark showed that the maximal Weissenberg number reached was no longer mesh-dependent [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF]. This idea was a new starting point and many improved numerical computations of viscoelastic flows were then performed (see e.g. [START_REF] Hao | Simulation for high Weissenberg number viscoelastic flow by a finite element method[END_REF][START_REF] Damanik | A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Afonso | The finite-volume method in computational rheology[END_REF] and [START_REF] Ciarlet | Handbook of numerical analysis[END_REF] for a recent review).

The main objective of the present paper is to bring some novelties in the challenging field of numerical methods for viscoelastic fluid flows. These novelties develop in three main axes:

1. A new and different log-conformation formulation for viscoelastic models is proposed. This formulation is non-singular when the Weissenberg number vanishes, while the original one, as proposed by Fattal and Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], degenerates, due to the apparition of the inverse of the Weissenberg number in the set of equations. With the present formulation, the problem reduces nicely to the Navier-Stokes equations at zero Weissenberg number. This is a major advantage, as it opens the door to continuation methods, starting smoothly at zero Weissenberg number and progressively increasing.

2. The steady problem is directly solved by a Newton method, while, to our knowledge, all previous approaches with the log-conformation formulation were time-dependent (see e.g. [START_REF] Damanik | A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems[END_REF]).

The derivation of a robust Newton solver for the strongly nonlinear steady log-formulation is based on an exact computation of all the derivatives. Notice that previous Newton methods relied on some finite difference methods for computing the Jacobian matrix, as the strong non-linearities were considered non-differentiable [START_REF] Kane | A comparison of four implementations of the logconformation formulation for viscoelastic fluid flows[END_REF][START_REF] Damanik | A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems[END_REF].

3. The discretization is based on an incompressible finite element method: the discrete velocity field satisfies the divergence free relation exactly. This is a major advantage when dealing with a transport equation: Recall that a non-zero divergence velocity field u introduces an additional source term div u in the stress transport equation: this term can then generate an exponential growth of the stresses, and then overrules all the benefits of using a logconformation formulation.

The outline of the paper is as follow. Section 1 starts with a presentation of the Johnson-Segalman model and its reformulation using the conformation tensor. Then the new non-singular logconformation formulation is presented. The first section ends with the variational formulation of the steady problem. Section 2 begins with the space discretization using an incompressible finite element method. Then, the upwind discretization of the transport term, based on the discontinuous Galerkin method is presented. Then, the Newton method applied to the discrete problem is introduced with special care to the derivation of highly nonlinear terms produced by the log-conformation formulation. Section 2 ends with the Euler-Newton continuation algorithm. Section 3 presents numerical computations on the bi-dimensional driven cavity flow. A specific adaptive mesh procedure allows resolving accurately both boundary layers and secondary vorties.

Finally, two appendices ellaborate on the computation of some strongly nonlinear terms and their derivatives.

1 Continuous setup

λ D a τ Dt + τ -2η p D(u) = 0 in ]0, t f [×Ω ρ Du Dt -div ( τ + 2η v D(u) -p.I) = 0 in ]0, t f [×Ω div u = 0 in ]0, t f [×Ω λ(τ -τ Γ ) = 0 on ]0, t f [×∂Ω - u = u Γ on ]0, t f [×∂Ω λτ (0) = λτ 0 and u(0) = u 0 in Ω where D Dt = ∂ ∂t + u.
∇u denotes the Lagrange derivative, λ 0 is the relaxation time, η p > 0 and η v 0 are the polymer and solvent viscosities, respectively, ρ 0 is the constant density, τ 0 , u 0 , τ Γ and u Γ are given initial and boundary conditions and ∂Ω -is the upstream boundary, defined by ∂Ω -= {x ∈ ∂Ω; u Γ .n(x) < 0}. The total Cauchy stress tensor is σ tot = -p I + 2η v D(u) + τ . When λ = 0, the problem reduces to the Navier-Stokes equations with a total viscosity η v + η p and τ = 2η p D(u).

Several studies published are concerned with the existence, uniqueness and regularity of a solution for this problem. In 1985, Renardy [START_REF] Renardy | Existence of slow steady flows of viscoelastic fluids with differential constitutive equations[END_REF] proved the existence of a solution to the steady problem. In 1987, Guillopé and Saut [START_REF] Guillopé | Résultats d'existence pour les fluides viscoélastiques à loi de comportement de type différentiel[END_REF] proved that there exists a unique strong solution local in time and that this solution is globally defined when the data are small and the fluid is not too elastic (see also [START_REF] Lin | On hydrodynamics of viscoelastic fluids[END_REF] for a different approach). In 1998, Fernández-Cara, Guillén and Ortega [START_REF] Fernández-Cara | Some theoretical results concerning non-Newtonian fluids of the Oldroyd kind[END_REF] extended this result to the large elasticity case and for an arbitrarily large but finite final time t f . In 2000, Lions and Masmoudi [START_REF] Lions | Global solutions for some Oldroyd models of non-Newtonian flows[END_REF] obtained an improved result (t f could be infinite) in the particular case of a corotational derivative (a = 0). In 2005, Kupferman, Mangoubi, and Titi [START_REF] Kupferman | A Beale-Kato-Majda breakdown criterion for an Oldroyd-B fluid in the creeping flow regime[END_REF] obtained a more explicit condition for the global in time existence of the solution when Ω = R 3 .

The solution (τ, u, p) exhibits some important properties. Hulsen [32, p. 6, eqn(30)] (see also Kwon and Leonov [38,p. 31,eqn (10)], Lee et al. [39, p. 383]) introduced the following definition of the conformation tensor when both a and λ are non-zero:

c = τ + η p aλ I
where I denotes the d × d identity matrix. Notice that c has the dimension of a stress and is undefined when λ = 0 or a = 0. Also, the conformation tensor c admits an explicit expression in integral form involving the unknown velocity field u [39, p. 385, eqn (3.29)]. From (1) we get

D a I Dt = -2aD(u) then D a τ Dt = D a c Dt - 2η p λ D(u).
Then, the Johnson-Segalman problem becomes:

(P 2 ): find c, u and p, defined in ]0, t f [×Ω, such that

λ D a c Dt + c = η p aλ I in ]0, t f [×Ω (3a) ρ Du Dt -div (c + 2η v D(u) -p.I) = 0 in ]0, t f [×Ω (3b) div u = 0 in ]0, t f [×Ω (3c) λ(c -c Γ ) = 0 on ]0, t f [×∂Ω - (3d) u = u Γ on ]0, t f [×∂Ω (3e) λc(0) = λc 0 and u(0) = u 0 in Ω (3f)
In 1990, Hulsen [START_REF] Hulsen | A sufficient condition for a positive definite configuration tensor in differential models[END_REF] showed that when the initial condition c 0 in (3f) is symmetric definite positive, and when ∇u is bounded, then c remains symmetric definite positive at any time t > 0.

Notice that when λ = 0, the problem degenerates, due to the apparition of aλ in the denominator in the constitutive equation (3a). In the case of the upper convected derivative (a = 1), Fattal and Kupferman [START_REF] Fattal | Constitutive laws for the matrix-logarithm of the conformation tensor[END_REF][START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF] introduced a change of variable ψ = log aλ η p c , the so-called log-conformation formulation. Nevertheless, when λ = 0 the expected solution is simply ψ = 0, while the reformulated problem still degenerates due to a division by λ in the conservation of momentum equation (see [17, p. 283], the last eqn of the page). This change of variable was then applied in [START_REF] Coronado | A simple method for simulating general viscoelastic fluid flows with an alternate log-conformation formulation[END_REF] to a generalized constitutive model covering the present Johnson-Segalman model (a ∈ [-1, 1]), but the obtained problem still degenerated when λ = 0. The next paragraph presents a slightly different change of variable that solves this degeneracy nicely. The new formulation will no longer be singular when λ = 0. This opens the possibility of building a robust steady solver based on the Newton method and a continuation algorithm, starting smoothly at λ = 0.

A new log-conformation formulation

In order to obtain a non-singular formulation, a different change of the logarithmic variable is introduced:

χ = η p aλ log aλ η p c = η p aλ log I + aλ η p τ (4a) ⇐⇒ c = η p aλ exp aλ η p χ and τ = η p aλ exp aλ η p χ -I (4b)
Notice that both χ, c and τ have the dimension of stress. At the limit case aλ = 0, the new variable χ is still well defined, as showed by a simple Taylor expansion of the matrix logarithm in the vicinity of aλ = 0:

χ = τ - aλ η p τ 2 + . . . + 1 n aλ η p n-1 τ n + . . .
For instance, when λ = 0, then τ = 2η p D(u) is the solution of the Johnson-Segalman model which reduces to a Newtonian fluid and χ = τ = 2η p D(u). For completeness, the derivation of the present log-conformation formulation is provided here, since there are some subtle modifications due to the different change of variable and the introduction of the the Gordon-Schowalter derivative parameter a. The impatient reader can jump directly to the new formulation (8a)-(8f).

The main technique of the present derivation of the log-conformation is related to the evolution of the principal axes of the conformation tensor [START_REF] Hulsen | Flow of viscoelastic fluids past a cylinder at high Weissenberg number: stabilized simulations using matrix logarithms[END_REF]. Recall that the symmetric positive definite matrix c can always be diagonalized as c = qcq T where c is diagonal and q is an orthogonal tensor, i.e. qq T = q T q = I. For convenience, let us denote by τ = Dτ Dt = ∂τ ∂t + u.∇τ the Lagrangian derivative of any tensor τ . Then ċ = q ċq T + qc qT + qcq T . Let us introduce r = q qT . By differentiating qq T = I we get qq T + q qT = 0 which also writes r T = -r. Then r is skew-symmetric.

We also obtain q = -rq and qT = -q T r T = q T r and then ċ = q ċq Tqcq T r Trqcq T . Next, let r = q T rq. Note that r is also skew-symmetric. We have r = qrq T and then

ċ = q ċ -cr T -rc q T (5a)
The generalized gradient introduced in (2) decomposes in this eigensystem as g a (u) = qg a q T where ga is not a priori diagonal, since in general c and g a (u) are not expected to share the same eigenvectors. Then -cg a (u) Tg a (u)c = q -cg T a -ga c q T and the Gordon-Schowalter derivative (1) writes:

D a c Dt = q ċ -c (r + ga ) T -(r + ga ) c q T (5b)
The constitutive equation (3a) leads to ċ = η p aλ 2 I -

1 λ c + (r + ga ) c + c (r + ga ) T (5c)
This is a tensorial equation. Let (c i ) 1 i d be the eigenvalues of c. Recall that c is diagonal and that its i-th diagonal entry is c i . Then, by taking the diagonal entries i = j, 1 i d, of the tensorial equation (5c), we get a differential equation for c i :

ċi = η p aλ 2 - 1 λ -2a di,i c i
where we have set d = q T D(u)q = ( di,j ) 1 i,j d . The previous relation also writes

ċ = η p aλ 2 I - I λ -2a diag( d) c (5d)
Now, let us perform the change of variable (4a)-(4b). The problem will be rewritten in terms of the variable

χ = η p aλ log aλ η p c = q η p aλ log aλ η p c q T = q χq T
where χ is the diagonal tensor with diagonal entries

χ i = η p aλ log aλ η p c i ⇐⇒ c i = η p aλ exp aλ η p χ i , 1 i d
Thus χi = η p aλ ċi ci and, since c, χ, ċ and χ are all diagonal:

χ = η p aλ ċc -1 = η p aλ η p aλ 2 I - I λ -2a diag( d) c c-1 , from (5d) = - 1 λ η p aλ I - aλ η p c-1 + 2η p λ diag( d)
From (4a), we have c -1 = aλ η p exp -aλ η p χ and then

χ = - 1 λ η p aλ I -exp - aλ η p χ + 2η p λ diag( d)
Then, replacing in (5a) the instances of c and c by χ and χ respectively, the Lagrangian derivative of χ reads:

χ = q χ -χr T -r χ q T = - 1 λ η p aλ I -exp - aλ η p χ + q 2η p λ diag( d) -χr T -r χ q T (5e)
Equivalently, this relation can be written

χ + φ a (χ, ∇u) + 1 λ η p aλ I -exp - aλ η p χ = 0 (5f)
where we have introduced the notation:

φ a (χ, ∇u) = q φa q T φa = χr T + r χ - 2η p λ diag( d)
Note that φ a (χ, ∇u) is symmetric since diag( d) is diagonal and r χ + χr T is symmetric. Notice also that both q and diag( d) are directly computable from χ and ∇u since q is the tensor containing the eigenvectors of χ and d = q T D(u)q. The computation of r χ + χr T is more technical. Let us compute r = (r i,j ) 1 i,j d . Since rT = -r we have ri,i = 0. By taking the off-diagonal entries i = j, 1 i, j d, in (5c), we get 0 = c i (-r i,j + ga;j,i ) + (r i,j + ga;i,j )c j where we have used rj,i = -r i,j and set ga;i,j = wi,j + a di,j and w = q T W (u)q = ( wi,j ) 1 i,j d . When c i = c j we obtain by simple development:

ri,j = ga;i,j c j + ga;j,i c i c i -c j = -wi,j + a c i + c j c i -c j di,j = -wi,j + a di,j tanh aλ η p χ i -χ j 2

Note that now r is directly computable from χ and ∇u but is still undetermined when χ i = χ j or aλ = 0. However, always when i = j and c i = c j :

φa;i,j = χ i rj,i + ri,j χ j = -(χ i -χ j ) ri,j = (χ i -χ j ) wi,j -     aλ η p χ i -χ j 2 tanh aλ η p χ i -χ j 2     2η p di,j λ As lim x→0 x tanh(x)
= 1, we obtain: φa;i,j = -2η p λ di,j when χ i = χ j or aλ = 0

Finally φ a (χ, ∇u) is directly computable from χ and ∇u, in all cases. Note that when all the eigenvalues of χ are equal, then φ a (χ, ∇u) = -2η p λ D(u). Also, when a = 0, we have

φ a (χ, ∇u) = χW (u) -W (u)χ - 2η p λ D(u).
This suggests the introduction of a nonlinear form κ(., .) such that

φ a (χ, ∇u) = χW (u) -W (u)χ + η p λ κ aλ η p χ, 2D(u) - 2η p λ D(u)
Then (5f) becomes

λ ( χ + χW (u) -W (u)χ) + η p κ aλ η p χ, 2D(u) + η p aλ I -exp - aλ η p χ = 2η p D(u) (5g) Let R d×d s denotes the set of d × d symmetric real matrix. The form κ(β, γ) is defined for all β, γ ∈ R d×d s by κ(β, γ) = qκq T (6a) κi,j = κ β i -β j 2 γi,j , 1 i, j d (6b)
where β = q T βq = diag(β i ) and γ = q T γq. Here, κ(x) is defined for all x ∈ R by

κ(x) = 0 when x = 0 1 - x tanh(x) otherwise (6c)
For convenience, let us define the following function, for all real µ 0 and

χ ∈ R d×d s f (µ, χ) =    0 when µ = 0 exp(µχ) -I µ -χ otherwise (7)
Note that the trace of λ f aλ η p , χ represents a local free energy stored by the micro-scale mechanical system (see e.g [START_REF] Bird | Dynamics of polymeric liquids[END_REF][START_REF] Wapperom | Thermodynamics of viscoelastic fluids: the temperature equation[END_REF][START_REF] Öttinger | Beyond equilibrium thermodynamics[END_REF]). This concept was recently used by Hu and Lelièvre [START_REF] Hu | New entropy estimates for the Oldroyd-B model and related models[END_REF][START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] for stability estimates and is employed later in this paper for an automatic adaptive meshing purpose.

The new log-conformation formulation of the Johnson-Segalman problem is:

(P 3 ): find χ, u and p, defined in ]0, t f [×Ω, such that

λ D 0 χ Dt + χ -f aλ η p , -χ + η p κ aλ η p χ, 2D(u) -2η p D(u) = 0 in ]0, t f [×Ω (8a) ρ Du Dt -div χ + f aλ η p , χ + 2η v D(u) -p.I = 0 in ]0, t f [×Ω (8b) div u = 0 in ]0, t f [×Ω (8c) λ(χ -χ Γ ) = 0 on ]0, t f [×∂Ω - (8d) u = u Γ on ]0, t f [×∂Ω (8e) λχ(0) = λχ 0 and u(0) = u 0 in Ω (8f)
where χ 0 , χ Γ , u 0 and u Γ are given. The functions f and κ are defined by ( 7) and (6a)-(6c), respectively. The elastic stress τ is computed explicitly as τ = χ + f aλ η p , χ and the total Cauchy stress tensor is

σ tot = -p I + 2η v D(u) + τ .
Remark 1 ( Limit λ = 0: Newtonian fluid) When λ = 0, the problem reduces as expected to the Navier-Stokes equations with a total viscosity η v + η p and χ = 2η p D(u).

Remark 2 ( Limit a = 0: corotational Johnson-Segalman fluid) When a = 0, the present log-formulation of the problem reduces nicely to the corotational Johnson-Segalman problem with τ = χ is the elastic stress, since f (0, .) and κ(0, .) both are zero, as showed in the next two properties.

Remark 3 ( Corotational versus Gordon-Schowalter tensor derivatives) Note the corotational derivative for χ in (8a). Thus, the present log-conformation formulation of the general Johnson-Segalman model can be interpreted as a nonlinear perturbation with nonlinear terms (f and κ) of the corotational Johnson-Segalman model. In the following section, we show that the numerical treatment of corotational derivative is much simpler than in the general Gordon-Schowalter case.

Property 1 ( Regularity of f ) The function f defined by ( 7) is continuously differentiable and

f (0, χ) = f (µ, 0) = 0, ∀µ ∈ R + , ∀χ ∈ R d×d s
Proof: This result follows from a simple Taylor expansion:

f (µ, χ) = µχ 2 I 2 + µχ 6 + . . . + (µχ) n (n + 2)! + . . . + Property 2 ( Skew symmetry of κ)
The function κ defined by (6a)-(6c) satisfies the following properties:

1. κ(., .) is nonlinear with respect to its first variable β and is linear with respect to its second variable γ.

2. κ is traceless:

tr κ(β, γ) = 0, ∀β, γ ∈ R d×d s
3. κ(β, γ) = 0 when β and γ are aligned, i.e. when they share the same eigensystem. As a special case, κ(β, γ) = 0 when all eigenvalues of β are equal.

4. κ is skew-symmetric with respect to its first variable:

κ(β, γ) : β = 0, ∀β, γ ∈ R d×d s
The skew symmetry extends to κ(β, γ) : σ = 0, ∀β, γ, σ ∈ R d×d s and β and σ share the same eigensystem (9)

5. κ is continuously differentiable everywhere and κ(0, γ) = 0 and ∂κ ∂β (0, γ) = 0, ∀γ ∈ R d×d s Proof: Remark that, from (6a)-(6c), we have κk,k = 0 for all k, 1 k d which leads to the traceless property. When β and γ are aligned then γ is diagonal and then κi,j = 0. The skew-symmetry result can be shown by the following development:

κ(β, γ) : β = (qκq T ) : (q β q T ) = d-1 i,j,k,l,m=0 q i,j κj,k q l,k q i,m β m q l,m = d-1 j,k,m=0 δ j,m δ k,m κj,k β m since d-1 i=0 q i,j q i,m = δ j,m and d-1 l=0 q l,k q l,m = δ k,m = d-1 k=1 κk,k β k = 0 since κk,k = 0
When β = 0, then all eigenvalues of β are equal and then κ(0, γ) (point 3). The proof of the differentiability of κ is postponed to appendix A.2, as it requires some technical developments.

Variational formulation of the steady problem

Here, the variational formulation of the steady version of the log-conformation formulation (8a)-(8f) of the Johnson-Segalman problem is considered. The inertia term (u.∇)u in the conservation of momentum (8b) is neglected. This is a common assumption in such flow simulations and the main difficulty is associated with the nonlinear terms related to the elasticity λ 0. Let us introduce three functional spaces:

T = L 2 (Ω) d×d s for square summable symmetric tensors, V (u Γ ) = {v ∈ H 1 (Ω) d ; v = u Γ
on ∂Ω} for velocities with square summable gradients which satisfies the boundary condition, and

Q = L 2 0 (Ω) = {q ∈ L 2 (Ω)
; Ω q dx = 0} for square summable pressures with zero average value. Next, let the following forms be defined for all χ, ξ ∈ T ,

u, v ∈ H 1 (Ω) d p, q ∈ L 2 (Ω): t(u; χ, ξ) = λ Ω ((u.∇)χ) : ξ dx + λ ∂Ω max(0, -u Γ .n) χ : ξ ds +λ Ω (χW (u) -W (u)χ) : ξ dx + η p Ω κ aλ η p χ, 2D(u) : ξ dx (10a) l(ξ) = λ ∂Ω max(0, -u Γ .n) χ Γ : ξ ds (10b) a 0 (χ, ξ) = Ω χ : ξ dx - Ω f aλ η p , -χ : ξ dx (10c) b 1 (ξ, u) = -2η p Ω ξ : D(u) dx (10d) b 2 (χ, v) = - Ω χ : D(v) dx - Ω f aλ η p , χ : D(v) dx (10e) c(u, v) = 2η v Ω D(u) : D(v) dx (10f) d(u, q) = - Ω q div u dx (10g)
The variational formulation writes

(F V ): find (χ, u, p) ∈ T × V (u Γ ) × Q such that t(u; χ, ξ) + a 0 (χ, ξ) + b 1 (ξ, u) = l(ξ), ∀ξ ∈ T (11a) b 2 (χ, v) -c(u, v) -d(v, p) = 0, ∀v ∈ V (0) (11b) -d(v, p) = 0, ∀q ∈ Q (11c)
The following properties are fundamental for solving the discrete problem and will admit a finitedimensional counterpart in the next section.

Property 3 ( Skew-symmetry of t) For all χ ∈ T and u ∈ V (0) such that div u = 0, we have

t(u; χ, χ) = 0 (12)
Proof: Integrating by part, we have

Ω ((u.∇)χ) : χ dx = - Ω ((u.∇)χ) : χ dx - Ω |χ| 2 div u dx + ∂Ω |χ| 2 u.n ds
and since u is divergence-free and vanishes on the boundary, the first term of the right-hand side of (10a), which gives the expression of t(u; χ, χ), is zero. Next, from the skew-symmetry of W (u) and the symmetry of χ we have (χW (u) -W (u)χ) : χ = 0. From property (9) we have κ aλ η p χ, 2D(u) : χ = 0 and hence get [START_REF] Crochet | Elastic effects in die entry flows[END_REF].

Property 4 ( Positivity of a 0 ) For all χ ∈ L 2 (Ω) d×d s , we have

a 0 (χ, χ) 0 ( 13
)
Proof: From (10c) and the definition (7) of f , we have

a 0 (χ, χ) = η p aλ Ω I -exp - aλ η p χ : χ dx
Note that σ : τ = tr(στ ) for all symmetric tensors σ and τ . Then, when aλ = 0:

a 0 (χ, χ) = η p aλ Ω tr I -exp - aλ η p χ χ dx
As χ and exp -aλ η p χ share the same eigensystem, if µ is an eigenvalue of χ then g(µ) is an eigenvalue of χ -exp -aλ η p χ χ where g(µ) = µ -exp(-aλµ/η p )µ. An easy inspection of the variation of g shows that g(µ) 0 for all µ ∈ R. We then obtain the result [START_REF] Crochet | On numerical die swell calculation[END_REF].

Remark 4 ( Corotational versus Gordon-Schowalter tensor derivatives (cont.)) Multiplying ( 1) by τ , integrating and rearranging we get

Ω D a τ Dt : τ dx = 1 2 d dt Ω |τ | 2 dx -2a Ω tr(τ D(u) τ ) dx
When a = 0, there is no way of determining the sign of the last term on the right-hand side of the previous relation. The corresponding term -a(τ D(u) + D(u)τ ), supported by a = 0 in the tensor derivative (1), acts as a source term in any tensorial transport problem involving the general Gordon-Schowalter derivative, such as the initial and conformation formulations of the Johnson-Segalman problem. This term is responsible for the observed exponential growth of the stress tensor τ and the failure of numerical methods. There is no hope of obtaining either a skewsymmetry or a positivity property for such problems, and this is especially true for the popular upper-convected derivative (a = 1). On the contrary, the log-conformation formulation involves a corotational tensor derivative of χ in (8a), leading to the skew-symmetry of t and the positivity of a 0 . Properties 3 and 4 are definitive advantages of the present log-conformation formulation over the initial and the conformation ones.

Note also that, when initially using a corotational derivative (a = 0), there is no need anymore for the log-conformation formulation. In that case, the new formulation presented here nicely reduces to the initial one with χ = τ while all the extra nonlinear terms disappear.

2 Discretization and numerical resolution

Space approximation

The main idea here is to replace the three spaces T , X, Q by some finite dimensional counterparts

T h ⊂ T , X h ⊂ X, Q h ⊂ Q in the variational formulation (11a)-(11c)
. When λ = 0 (a Newtonian fluid), a finite dimensional linear system is obtained, whose matrix has the following bloc structure:

  A B T 0 B -C -D T 0 -D 0  
This linear system is often called the three field Stokes problem and the approximation space pairs (T h , X h ) and (X h , Q h ) both need to satisfy a compatibility condition, known as the inf-sup or Babuška-Brezzi condition [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[END_REF][START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. There are many possible choices for T h , X h and Q h , each of them having some advantages and drawbacks. In this paper, we consider the following choice (see macro-element technic P 2 continuous P 1 discontinuous u σ and p

T h = {τ h ∈ T ; τ h|K ∈ P 1 , ∀K ∈ T h } X h = {v h ∈ (H 1 (Ω) ∩ C 0 (Ω)) d ; v h|K ∈ P 2 , ∀K ∈ T h } V h (g) = X h ∩ V (π h (g)) Q h = {q h ∈ Q; q h|K ∈ P 1 , ∀K ∈ T h }
where T h is a finite element mesh of the flow domain Ω and h > 0 denotes the largest edge length of the mesh. Notice that both T h and Q h contain discontinuous piecewise polynomials while X h elements are continuous. The discontinuous approximation of stresses has two advantages. On one hand, D(X h ) ⊂ T h and then the compatibility condition between T h and X h is satisfied for the bi-linear form b 1 , as showed by Fortin and Pierre [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF][START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF]. On the other hand, it allows an efficient treatment of the stress transport term by the discontinuous Galerkin method.

The discontinuous approximation of pressure also presents a major advantage. As div(X h ) ⊂ Q h , it leads to an exact divergence-free approximation of the velocity: for any field v h ∈ X h satisfying d(v h , q h ) = 0 for all q h ∈ Q h , we have div v h = 0 point-wise everywhere in Ω. The pair (X h , Q h ) is known as the Scott-Vogelius lowest-order finite element approximation [START_REF] Scott | Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials[END_REF]. This is a major advantage when dealing with a transport equation. The only drawback is that the pair (X h , Q h ) does not satisfy the inf-sup condition for an arbitrary mesh. There exists a solution to this however: Arnold and Qin [START_REF] Arnold | Quadratic velocity/linear pressure Stokes elements[END_REF] proposed a macro element technique applied to the mesh [START_REF] Arnold | Quadratic velocity/linear pressure Stokes elements[END_REF] that allows satisfying the inf-sup condition: for any triangular finite element mesh, it is sufficient to split each triangle in three elements from its barycenter (see Fig. 1). This technique has been implemented for the present computations. Notice that the macro element technique extends to quadrilateral meshes [START_REF] Arnold | Quadratic velocity/linear pressure Stokes elements[END_REF] and to the three-dimensional case [START_REF] Zhang | A new family of stable mixed finite elements for the 3d Stokes equations[END_REF].

Let us review and discuss some other possible choices for the discretization of the problem.

• A very popular choice consists in using staggered finite difference grids for the velocitypressure pair. With this choice the approximate velocity fields are also exactly divergencefree, the finite-difference implementation is simple and this approximation also extends in the mixed finite element context with the incompressible Raviart-Thomas finite element [24, chap. 3]. For the stress approximation, there are several possibilities: in [57, p. 17], the present author also used a staggered grid approximation (diagonal stress components at cell centers and off-diagonal components at cell corners in the bi-dimensional case). With this choice, the stress-velocity pair satisfies the inf-sup condition [START_REF] Saramito | Numerical simulation of viscoelastic fluid flows using incompressible finite element method and a θ-method[END_REF] and leads to a robust scheme that is able to reach solutions at high Weissenberg numbers (see [START_REF] Saramito | Efficient simulation of nonlinear viscoelastic fluid flows[END_REF][START_REF] Saramito | Flow characteristics of viscoelastic fluids in an abrupt contraction by using numerical modeling[END_REF]). In [18, p. 27],

Fattal and Kupferman used a cell-centered approximation for all the stress components. As they pointed out [18, p.29]:

"This implies that the rotational components of the system may be sensitive to numerical instabilities. A natural remedy would have been to store the off-diagonal elements of the stress tensor at cell corners. This cannot be done in our framework as the log-conformation tensor is an entity whose tensorial nature is essential."

See also point 6 of the discussion, p. 36 of the same paper. With this choice, it is not clear whether the stress-velocity pair satisfies the inf-sup condition or not. Nevertheless, using this approximation, Fattal and Kupferman were able to compute solutions for some high Weissenberg numbers. A possible explanation for this success can be found in a paper by Baranger and Sandri [START_REF] Baranger | Finite element approximation of viscoelastic fluid flow: existence of approximate solutions and error bounds. part i. discontinuous constraints[END_REF], who have shown that the stress-velocity compatibility condition exhibited by [START_REF] Fortin | On the convergence of the mixed method of Crochet and Marchal for viscoelastic flows[END_REF] is only required when η v /(η v + η p ) ≪ 1, i.e. in the absence of the pure viscous contribution. Otherwise, when η v is not too small with respect to η p , a much larger class of discretization schemes is allowed. This last condition was fulfilled by Fattal and Kupferman: all the computations presented by these authors used η v = η p .

• A classical choice is to use the Taylor-Hood P 2 -P 1 continuous element [START_REF] Taylor | A numerical solution of the Navier-Stokes equations using the finite element technique[END_REF] (see also [9, p. 252]) for the velocity-pressure spaces (X h , Q h ) together with piecewise linear and discontinuous stresses. This choice of stress approximation was first introduced by A. and M. Fortin [START_REF] Fortin | A new approach for the FEM simulation of viscoelastic flows[END_REF]. Damanik [14, p. 25] recently used a variant for quadrilaterals meshes, combining Q 2 continuous velocities with P 1 discontinuous pressures [9, p. 216] and a quadratic continuous approximation for the stresses. As div(X h ) ⊂ Q h , for these elements, the approximate velocity field are not exactly divergence-free. This is the main difference with the Scott-Volelius element used in the present computation. We performed tests for both the Taylor-Hood approximation and the Q 2 -P 1,d variant for quadrilaterals, together with P 1,d stresses, on the driven cavity problem and observed that the steady approximate solution presents some problems: for large dimensionless Weissenberg numbers (e.g. W e = 1), the divergence div u h does not converge anymore to zero with mesh refinement (h → 0). This could be due to the low regularity of the velocity field in the boundary layers of the cavity, especially at high Weissenberg numbers. This motivates the use the divergence-free approximation for the velocity-pressure pair.

• Another possible choice is to use constant and discontinuous P 0 stresses together with the divergence-free Scott-Vogelius velocity-pressure space pair. Notice that we have now D(Xh) ⊂ T h. This choice was first suggested by Mangoubi and Boyaval [START_REF] Mangoubi | Analytical and numerical problems in the computation of the flow of viscoelastic fluids[END_REF][START_REF]Modélisation mathématique et simulation numérique en science des matériaux[END_REF][START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF] in a theoretical paper, who were able to prove a nice stability property for the free energy. There were no numerical experiment with this combination. Moreover, we also performed some test on the driven cavity problem and observed that this choice is less robust and less precise than choosing piecewise linear and discontinuous stresses as done here.

This discretization review is far from complete and exhaustive: there are many other possible approaches and some of them have already been tested. Nevertheless, note that the exact divergencefree property for the approximation of the velocity field appears as an essential condition in such flow simulations. Both the staggered finite difference method and the Scott-Vogelius finite element one for the velocity-pressure pair satisfy this condition. The Scott-Vogelius finite element method presents two additional advantages compared to the finite difference one: (i) it allows the same approximation for all the stress components while satisfying the inf-sup condition and (ii) it is more flexible when dealing with complex geometries or adaptive meshes.

Approximation of the transport term

Let us turn to the discretization of the nonlinear stress transport term by the discontinuous Galerkin method. We introduce the following form, defined for all u ∈ H 1 (Ω) d and χ, ξ ∈

L 2 (Ω) d×d s such that χ |K , ξ |K ∈ H 1 (K) d×d s for all element K ∈ T h : t h (u; χ, ξ) = λ K∈T h K ((u.∇)χ) : ξ dx + λ S∈S (i) h S [[χ]] : θ 2 |u.n| [[ξ]] -(u.n) { {ξ} } ds + λ ∂Ω max (0, -u Γ .n) χ : ξ ds + λ Ω (χW (u) -W (u)χ) : ξ dx + Ω κ aλ η p χ, 2D(u) : ξ dx ( 14
)
where

S (i) h
denotes the set of internal sides of the mesh T h , [[ξ]] is the jump of a piecewise discontinuous tensor across a side and { {ξ} } is its average value [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF][START_REF] Saramito | Efficient C++ finite element computing with Rheolef: volume 2: discontinuous Galerkin methods[END_REF]. The three first terms represent a discrete counterpart of the transport term that can not be defined globally since the stress approximation is piecewise discontinuous. A term is weighted by a factor θ 0. Choosing θ = 0 corresponds to the so-called centered flux approximation, while θ > 0 is the upwind flux approximation. The case of θ = 1, the most popular upwind discontinuous approximation scheme, is considered here. The upwind technique is known as an efficient approach to avoid spurious oscillations in the approximate solutions.

The discrete variational formulation of the problem writes:

(F V ) h : find χ h ∈ T h , u h ∈ V h (u Γ ) and p h ∈ Q h such that t h (u h ; χ h , ξ h ) + a 0 (χ h , ξ h ) + b 1 (ξ h , u h ) = l(ξ h ), ∀ξ h ∈ T h b 2 (χ h , v h ) -c(u h , v h ) -d(v h , p h ) = 0, ∀v h ∈ V h (0) -d(u h , q h ) = 0, ∀q h ∈ Q h
where the continuous tri-linear form t has simply been replaced by t h and the functional spaces by their finite dimensional counterpart.

Property 5 ( Discrete generalized skew-symmetry) For all χ h ∈ T h and u h ∈ V h (0) such that div u h = 0, we have

t h (u h ; χ h , χ h ) 0 ( 15
)
with an equality when θ = 0.

Proof: Integrating by part on an element K leads to

K ((u h .∇)χ h ) : χ h dx = - K χ h : ((u h .∇)χ h ) dx + K |χ h | 2 (u h .n) dx since div u h = 0. Then K∈T h K ((u h .∇)χ h ) : χ h dx = 1 2 K∈T h ∂K |χ h | 2 (u h .n) ds = 1 2 S∈S (i) h S [[|χ h | 2 ]] (u h .n) ds
Note that for any discontinuous scalar field φ and ϕ across a side S we have

[[φϕ]] = [[φ]]{ {ψ} } + { {φ} }[[ϕ]]. Thus [[|χ h | 2 ]]/2 = [[χ h ]] : { {χ h } } and K∈T h K ((u h .∇)χ h ) : χ h dx = S∈S (i) h S [[χ h ]] : { {χ h } } (u h .n) ds
Using u h = 0 on ∂Ω and dealing with the two last terms of ( 14) as in the proof of property 3, we obtain:

t h (u h ; χ h , χ h ) = θλ 2 S∈S (i) h S |[[χ h ]]| 2 |u h .n| ds
which completes the proof.

Property 6 ( Discrete positivity) For all χ h ∈ T h , we have a 0 (χ h , χ h ) 0 ( 16)

Proof: As T h ⊂ L 2 (Ω) d×d s
, this result is a direct consequence of property 4.

Newton method

Newton methods for the numerical resolution of steady viscoelastic flow problems started in the early eighties with the work of Bézy, Crochet and Keunings [START_REF] Crochet | Elastic effects in die entry flows[END_REF][START_REF] Crochet | On numerical die swell calculation[END_REF]. In 1992, Fortin and Zine [START_REF] Fortin | An improved GMRES method for solving viscoelastic fluid flow problems[END_REF] proposed a quasi-Newton variant, where the Jacobian matrix was approximated instead of being completely recomputed at each iteration. After a long time, the Newton method approach for viscoelastic fluid problems is now coming back. While previous works relied on some finite difference methods for computing the Jacobian matrix, in 2009, Howell [START_REF] Howell | Dual-mixed finite element approximation of Stokes and nonlinear Stokes problems using trace-free velocity gradients[END_REF] computed exactly the linearized steady Johnson-Segalman problem in its initial formulation. In 2010, Damanik et al [START_REF] Damanik | A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems[END_REF][START_REF] Damanik | FEM simulation of non-isothermal viscoelastic fluids[END_REF] turned to a Newton method for the time-dependent log-conformation formulation. They used a finite difference method for computing the Jacobian matrix. Let us quote a remark made in the 2009 paper by Kane, Guénette, and Fortin [35, p. 45] for the log-conformation formulation: "[...] there is no hope to fully linearize the constitutive equations for the Newtons method without using some numerical tricks such as finite difference methods. The bulk of the computations relies indeed on the calculation, at each Gauss node, of the eigenvalues and eigenvectors of the conformation tensor which are not differentiable functions".

In the present paragraph, an exact expression of the derivatives for the fully linearized constitutive equation is presented for the log-conformation formulation, without any trick such as finite difference methods for computing the Jacobian matrix. The present Newton method directly treats the steady problem: by an obvious extension, it applies also to fully implicit time dependent simulations.

The discrete problem can be put in a compact form:

find (χ h , u h , p h ) ∈ T h × V h (u Γ ) × Q h such that F (λ; (χ h , u h , p h )) = 0
where F is defined in variational form for all (ξ h , v h , q h ) by

F (λ; (χ h , u h , p h )), (ξ, v, q) = t h (u h ; χ h , ξ h ) + a 0 (χ h , ξ h ) + b 1 (ξ h , u h ) -l(ξ h ) + b 2 (χ h , v h ) -c(u h , v h ) -d(v h , p h ) -d(u h , q h )
and where ., . stands for the duality product induced by the L 2 pivot space, i.e. ϕ, φ = Ω ϕφ dx for all ϕ, φ defined in Ω. The function F has two variables λ ∈ R + and

U = (χ h , u h , p h ) ∈ T h × X h × Q h .
The λ variable will be used as a continuation parameter in the next section. The Newton method defines the sequence χ

(k) h , u (k) h , p (k) k 0
by recurrence as:

• k = 0: let χ (0) h , u (0) h , p (0) ∈ T h × V h (u Γ ) × Q h being given. • k 0: let χ (k-1) h , u (k-1) h , p (k-1) ∈ T h × V h (u Γ ) × Q h being known. Find (δχ h , δu h , δp h ) ∈ T h × V h (0) × Q h such that ∂F ∂U λ; χ (k-1) h , u (k-1) h , p (k-1) h .(δχ h , δu h , δp h ) = -F λ, χ (k-1) h , u (k-1) h , p (k-1) h
and then defines

χ (k) h = χ (k-1) h + δχ h , u (k) h = u (k-1) h + δu h and p (k) h = p (k-1) h + δp h
At each step k 0, this algorithm solves a linear subproblem involving the Jacobian ∂F ∂U . The Newton method has only local convergence properties, i.e. the initial value should be close enough to the solution. In order to circumvent this limitation, a globalized Newton variant is used here. It is based on a damping strategy, as described and implemented in the Rheolef free software FEM library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF]. Note the absolute value function that appears in the upwind term in the definition [START_REF] Damanik | FEM simulation of non-isothermal viscoelastic fluids[END_REF] of t h :

s h (u h ; χ h , ξ h ) = 1 2 S∈S (i) h S [[χ h ]] : [[ξ h ]] |u h .n| ds
This term is not differentiable with respect to u h and neither is t h nor F . Nevertheless, the absolute value is convex and some tools from subdifferential calculus can be used to overcome this new difficulty. Let us introduce the multi-valued sign function:

sgn(x) =    {1} when x > 0 [-1, 1] when x = 0 {-1} when x < 0
Then, the subdifferential of the absolute value function is sgn(x) and for all δw h , w h , u h , v h ∈ X h , we introduce a generalization of the partial derivative:

∂s h ∂u h (u h ; χ h , ξ h ).(δu h ) = 1 2 S∈S (i) h S [[χ h ]] : [[ξ h ]] sgn(u h .n) δu h .

n ds

The Jacobian ∂F ∂U can thus be defined as a multi-valued subgradient set. In 1993, Qi and Sun [START_REF] Qi | A nonsmooth version of Newton's method[END_REF] showed how the Newton method extends to this case: any element that belongs to this set represents a valid direction for the Newton correction step. Such a non-smooth Newton method was successfully implemented for the steady Navier-Stokes equations with the discontinuous Galerkin method and upwind scheme [START_REF] Saramito | Efficient C++ finite element computing with Rheolef: volume 2: discontinuous Galerkin methods[END_REF].

The multi-valued Jacobian

∂F ∂U is defined for all (χ h , u h , p h ) ∈ T h × V (u Γ ) × Q h , and
(δχ h , δu h , δp h ) ∈ T h × V (0) × Q h (ξ h , v h , q h ) by ∂F ∂U (λ, (χ h , u h , p h )) .(δχ h , δu h , δp h ), (ξ h , v h , q h ) = a 1 (χ h , u h ; δχ h , ξ h ) + b 11 (χ h , u h ; ξ h , δu h ) + b 12 (χ h ; δχ h , v h ) -c(δu h , v h ) -d(v h , δp h ) -d(δu, q)
where a

+ Ω κ aλ η p χ, 2D(u) : ξ dx b 12 (χ; δχ, v) = ∂b 2 ∂χ (χ, v).(δχ) = - Ω δχ : D(v) dx - Ω exp ′ aλ η p χ : δχ : ξ dx
and where t 0,h represents

t 0,h (u h ; χ h , ξ h ) = K∈T h K ((u h .∇)χ h ) : ξ h dx - S∈S (i) h S [[χ h ]] : { {ξ h } } (u h .n) ds + Ω (χ h W (u h ) -W (u h )χ h ) : ξ h dx
Here exp ′ (χ) : δχ denotes the differential at δχ of the exponential of a d × d matrix. Recall that when χ and δχ commute, then exp ′ (χ) : δχ = exp(χ) δχ while the general case is more complex. Also, ∂κ ∂β involves the derivatives of eigenvalues and eigenvectors of a d × d matrix with respect to the matrix coefficients. Note that since κ is linear with respect to its second variable γ we have ∂κ ∂γ (β, γ) : δγ = κ(β, δγ). A key point of the present implementation of the Newton method is the exact computation of both exp(χ), exp ′ (χ), κ(β, γ) and ∂κ ∂β (β, γ). Indeed, these computations can be performed explicitly: Rouvière [56, p. 297] proposes some tools for the derivation of the exponential of a matrix, while Hairer and Wanner [27, p. 102] give some formula for computing the derivative of eigenvalues and Magnus and Neudecker [43, p. 177] also present the derivative of eigenvectors with respect to the coefficient of a matrix. Kane, Guénette, and Fortin [35, p. 50] presented some explicit expressions for exp(χ) and exp ′ (χ) based on the mapple software: these expressions contain some errors and we provide the correct expressions in appendix A.1 for completeness. Appendix A.2 shows that κ is continuous and differentiable with respect to the matrix coefficients and details for the first time the computation of these quantities. Finally, the integrals involving nonlinear expressions are evaluated by a Gauss quadrature formula with six interior nodes in a triangle: this formula is exact for polynomials whith degree lower or equal to four. Numerical experiments with higher order of quadrature do not show significant changes in the numerical results.

Euler-Newton continuation algorithm

The aim of the Euler-Newton continuation algorithm is to start from a previously computed solution at a smaller λ (or dimensionless Weissenberg number), perform a prediction using an Euler scheme (as if λ was a pseudo-time) and then perform corrections using the Newton method. This approach allows to reach high values of the λ parameter efficiently.

Let us denote U = (χ h , u h , p h ) ∈ T h × X h × Q h , such that the approximate nonlinear problem can be written in a concise form F (λ, U ) = 0. In this paragraph, we introduce an Euler-Newton continuation algorithm (see e.g. [52, p. 176] or [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]) that consists of two nested loops: algorithm 1 (continuation)

• n = 0: Let (λ 0 , U 0 ) be given. Compute U0 = - ∂F ∂U (λ 0 , U 0 ) -1 ∂F ∂λ (λ 0 , U 0 )
• n 0: Assume (λ n , U n ) and Un are known.

1) First choose a step ∆λ n and set λ n+1 = λ n + ∆λ n .

2) Then, perform an Euler prediction by computing

w 0 = U n -∆λ n ∂F ∂U (λ n , U n ) -1 ∂F ∂λ (λ n , U n )
3) Then, perform a Newton correction step: for all k 0, with W k being known, compute

W k+1 = W k - ∂F ∂U (λ n+1 , W k ) -1 F (λ n+1 , W k )
Once the the correction loop has converged, set U n+1 = W ∞ .

4) Finally, compute

Un+1 = - ∂F ∂U (λ n+1 , U n+1 ) -1 ∂F ∂λ (λ n+1 , U n+1 )
The step ∆λ n can either be fixed or chosen by adjusting the contraction ratio of the Newton method [START_REF] Paumier | Bifurcation et méthodes numériques. Applications aux problèmes elliptiques semi-linéaires[END_REF][START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]. The previous algorithm requires the computation of ∂F ∂λ , which is given by:

∂F ∂λ (λ, (χ h , u h , p h )) , (ξ h , v h , q h ) = t 0,h (u h ; χ h , ξ h ) + s h (u h ; χ h , ξ h ) + ∂Ω max(0, -u Γ .n) χ h : ξ h ds + a η p Ω ∂κ ∂χ aλ η p χ h , 2D(u h ) : χ h : ξ h dx - a η p Ω ∂f ∂µ aλ η p , -χ h : ξ h dx - a η p Ω ∂f ∂µ aλ η p , χ h : D(v h ) dx where ∂f ∂µ (µ, χ) =          χ 2 2 when µ = 0 - 1 µ 2 (exp (µχ) -I) + 1 µ exp ′ (µχ) : χ otherwise 2.

Automatic adaptive mesh

The anisotropic auto-adaptive mesh feature available in the Rheolef free software FEM library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF] has been used in this paper. This feature is based on the free software bamg bidimensional anisotropic mesh generator developed by Hecht [START_REF] Hecht | BAMG: bidimensional anisotropic mesh generator[END_REF]. Let us first summarize the principle of the adaptive mesh procedure. Let T h be an initial mesh and U = (χ h , u h , p h ) ∈ T h × X h × Q h be the solution of the discrete nonlinear problem (F V ) h associated to the mesh T h . Let φ be a governing field to be suitably chosen from the solution U . For a piecewise linear interpolation of φ, the interpolation error in the unitary direction d ∈ IR 2 is estimated at any element K ∈ T h by:

e K,d = h 2 K,d ∂ 2 φ ∂d 2 in K
where h K,d is the length of K in the d direction and

∂ 2 φ ∂d 2 = d T H(φ)d and H(φ) =     ∂ 2 φ ∂x 2 ∂ 2 φ ∂x∂y ∂ 2 φ ∂x∂y ∂ 2 φ ∂y 2    
Here H(φ) denotes the Hessian of φ. One possible way of adapting the mesh to the computation of φ is to equi-distribute this error, i.e. to make it constant over all triangles and in all directions. Let λ 1 , λ 2 be the eigenvalues of H(φ) and d 1 and d 2 the associated eigenvectors:

∂ 2 φ ∂d 2 1 = λ 1 and ∂ 2 φ ∂d 2 2 = λ 2
The error e K,d is independent of d and K when there exists a constant e 0 > 0 independent of K such that e K,d1 = e K,d2 = e 0 . This writes equivalently

F e 1 e 2 h 1 d 1 d 2 h 2
h 2 K,d1 |λ 1 | = h 2 K,d2 |λ 2 | = e 0 , ∀K ∈ T h
The constant e 0 represents a global surface density factor. The adapted mesh generator tries to shrink elements in all directions with a factor √ e 0 , hence all the bidimensional element areas are reduced by a factor e 0 . Suppose that H(φ) is non-singular, i.e. λ 1 λ 2 = 0. The constant e 0 being known, our aim is now to build triangles of length h i in the d i direction with h i = e 0 /|λ i |, i = 1, 2. Such a triangle has no privileged direction in a metric such that the two h i d i vectors, i = 1, 2, have the same norm. Thus, let us introduce the metric M(φ) tensor that have the same eigenvectors as H(φ) and |λ i |, i = 1, 2 as eigenvalues. The induced norm . M satisfies

h i d i M = h K d T i M(φ)d i = √ e 0 , i = 1, 2
It then suffices to build an isotropic mesh in the Riemann space associated to the metric M(φ).

In the Euclidean space, this mesh locally shrinks with a factor h i in the d i direction.

It remains to choose a suitable φ governing field. Several numerical experiments have motivated the use of a combination of a free energy and of a viscous dissipation term:

φ = λ tr f aλ η p , χ + η v |D(u)| 2
where f is defined in [START_REF] Boyaval | Free-energy-dissipative schemes for the Oldroyd-B model[END_REF]. As both χ h and D(u h ) are linear and piecewise discontinuous, the governing field φ is also approximated by a piecewise linear and discontinuous function φ h . Obtaining its hessian thus requires a L 2 projection and two derivations. The L 2 projection of φ h gives a piecewise linear and continuous approximation of φ, denoted by φh . The discrete Hessian H h is obtained from φh by computing first its discrete gradient g h , continuous and piecewise linear, from the variational formula

Ω g h .v h dx = Ω ∇ φh .v h dx
where v h is any continuous piecewise linear vector. Then H h is obtained from the variational formula

Ω H h : ξ h dx = Ω D(g h ) : ξ h dx
where ξ h is any continuous piecewise linear tensor. Solving a problem using such an automatic adaptive mesh is an iterative process, which involves three main steps :

1. Starting from an initial mesh T h , solve the problem using the Newton method. Let U be the corresponding solution associated to the mesh T h .

2. From U , compute the governing field φ h .

3. From the governing field φ h , defined on the mesh T h , generate a totally new mesh, denoted by T

Then, T

is used to solve the problem, and so on, until convergence of both the mesh and its associated solution. Hence, the final adapted mesh minimizes the interpolation error for the governing field. This choice of the governing field has been found able to catch accurately the boundary layers and secondary vorticies, as showed in the next section.

Tests on the smoothed driven cavity flow

For the purpose of comparison with previous studies [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF][START_REF] Hao | Simulation for high Weissenberg number viscoelastic flow by a finite element method[END_REF], we consider the steady bi-dimensional smoothed driven cavity benchmark problem with Ω = ]0, L[ 2 with L > 0. The boundary velocity u Γ is zero except on the top boundary {y = L} where u Γ (x, L) = (16U x 2 (L -x) 2 /L 2 , 0) with U > 0. The fluid parameters are a = 1 and η p = η v and the dimensionless Weissenberg number is W e = λU/L. Despite its simple geometry, this is a very difficult problem: to our knowledge, most numerical methods based on the initial formulation failed for W e 0.1. Due to singularities near corners between the lid and the sidewalls, the lid-driven cavity may encounter the high W e number problem or present many challenges in the form of singularity points in the flow. Therefore, it is known as a highly stringent test problem for numerical methods. Fig. 3.a plots the residual term F (λ, U ) L 2 versus the iteration number n iter of the damped Newton method. The loop is initialized from the Newtonian solution, associated to W e = 0 and the solution is then computed directly for a specific W e. Observe the quadratic convergence in log scale, up to W e = 0.5. The algorithm stops when the machine precision is reached. For the last test case, W e = 0.6 and the damped strategy still allows the convergence of the algorithm at the price of a larger number of iteration and the loss of the quadratic convergence. In that case, the Euler-Newton continuation algorithm, which restarts from a previous computed solution at a smaller W e, is more efficient. We used a Weissenberg number increment of ∆W e = 0.1 in the present computations. Moreover, the Euler-Newton continuation algorithm is more robust and allows reaching solutions at much higher Weissenberg numbers. With this algorithm, each step uses about five resolutions of the Jacobian linear system and the solution at W e = 1 can be reached from the solution at W e = 0 with about fifty resolutions of linear systems. This convergence property is mesh-invariant, as showed by Fig. 3.b (for the mesh-invariance property of nonlinear algorithms, see [61, chap. 8]). This algorithm provides also all the intermediate solutions at W e = 0.1, . . . , 0.9. Such very efficient viscoelastic computations for large Weissenberg numbers can be compared with the thousands of steps involving linear systems and required by all the current time-dependent approaches for reaching the steady solution [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF][START_REF] Hao | Simulation for high Weissenberg number viscoelastic flow by a finite element method[END_REF][START_REF] Kane | A comparison of four implementations of the logconformation formulation for viscoelastic fluid flows[END_REF][START_REF] Damanik | A monolithic FEM approach for the logconformation reformulation (LCR) of viscoelastic flow problems[END_REF]. Moreover, the number of time steps required by these methods increases with the Weissenberg number, and also with the mesh refinement, when non-fully implicit time schemes are used.

The current implementation uses the Rheolef free software FEM library [START_REF] Saramito | Efficient C++ finite element computing with Rheolef[END_REF][START_REF] Saramito | Efficient C++ finite element computing with Rheolef: volume 2: discontinuous Galerkin methods[END_REF] that is available as a standard package under the Debian and Ubuntu GNU/Linux systems. The Jacobian matrix ∂F ∂U is large and sparse: the linear system is solved with a direct method with the help of the parallel and massively distributed memory linear solver mumps [START_REF] Amestoy | Multifrontal massively parallel solver (MUMPS). User's guide[END_REF] together with the scotch [START_REF] Pellegrini | PT-Scotch and libscotch 5.1 user's guide[END_REF] ordering algorithm for minimizing the fill-in of the sparse matrix. For the largest meshes (h = 1/80, 1/160 and the adaptive meshes), the computations are run on 32 processors on a BullX DLC supercomputer (Bull Newsca) composed of nodes having two Intel Sandy-Bridge processors and connected to a FDR infiniband non-blocking low latency network. The computation of a whole branch of solutions is performed in less than one hour of real time.

Table 1 (left) groups the maximal Weissenberg number, denoted as W e max , reached by the Euler-Newton continuation algorithm: for higher W e, the algorithm stops, due to a singular solution where the Jacobian matrix ∂F ∂U is non-invertible. When such a singularity occurs in the continuation algorithm, the Weissenberg number increment is divided by two and the iteration restarted. The Weissenberg number increment is limited to 10 -7 . When using an adaptive mesh, we observe that the Euler-Newton continuation algorithm is able to jump behind the singular point and to continue its progression on a branch of solutions until a second singular point is reached near h W e max 1/20 W e = 3. In that case, the determinant of the Jacobian matrix exhibits a change of sign after the first singular point (see e.g. [START_REF] Seydel | Practical bifurcation and stability analysis[END_REF]). The Weissenberg number where this change of sign occurs is denoted as W e c on Table 1 (right). Note that for the finest uniform mesh (h = 1/160), the singular point is W e c = 2.08 and for the finest adapted one, W e c = 2.03. Then W e c seems to become mesh insensitive with mesh refinement. A deeper analysis of singular points is planed as future work. It requires more advanced tools than the simple continuation Euler-Newton algorithm. The first and second singular points could be e.g. bifurcation points associated to a loss of stationarity of the solution: for larger W e, the solution becomes non-stationary, as suggested by time-dependent simulations (see e.g. [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF], Fig. 8). Pakdel, Spiegelberg and McKinley also showed by experimental observations [START_REF] Pakdel | Cavity flows of elastic liquids: twodimensional flows[END_REF] for the motion of viscoelastic fluids in the lid-driven cavity geometry that, at large Weissenberg numbers, the fluid motion becomes unstable and a three-dimensional flow develops. Figs. 5 and6 show the adapted meshes and stream isovalues of the function for W e = 1, 2 and 3. The stream function ϕ is defined as the unique function that satisfies -∆ϕ = ∂ y u x -∂ x u y in Ω with ϕ = 0 on ∂Ω. Ten negative and fifteen positive equi-spaced isolines are represented on each plot. Observe that viscoelastic effects break the symmetry observed for the velocity field of cavity flows of viscous Newtonian fluids at zero Reynolds number. At low Weissenberg number, the flow remains two-dimensional but the center of the primary recirculating vortex in the cavity shifts progressively upstream (left). These results are in qualitative agreement with experimental results [START_REF] Pakdel | Cavity flows of elastic liquids: twodimensional flows[END_REF]. Note also that inertia in viscous Newtonian fluids (i.e. W e = 0 and Re > 0) has opposite effects (see e.g. [62, part 2]): the center of the primary recirculating vortex in the cavity shifts downstream (right).

x Table 3: The Oldroyd-B problem on the driven cavity with η p = η v . Dimensionless main and secondary vortices center position and activity for W e = 1, 2 and 3.

Table 3 groups the main and secondary vortices center position and activity for W e = 1, 2 and 3. These values are provided for future cross validation purpose. One can observe that the main vortex goes left and decreases in activity while the center of the left secondary vortex goes up and right and its activity increases. Also, the center of the right secondary vortex roughly stays in place while its activity remains constant. Finally, Fig. 9 plots the cut of the first component u x of the velocity and the χ xx component along the x = 1/2 vertical axis, as computed with the adaptive mesh. The computations with the finest uniform mesh (h = 1/160) give similar results (not shown as the difference is not graphically perceptible). Observe the excellent quantitative agreement with both the computations of Hao and Pan [START_REF] Hao | Simulation for high Weissenberg number viscoelastic flow by a finite element method[END_REF] and Fattal and Kupferman [START_REF] Fattal | Time-dependent simulation of viscoelastic flows at high Weissenberg number using the log-conformation representation[END_REF].

Conclusion

The new log-conformation formulation of viscoelastic fluid flow presented in this paper allows a direct steady numerical resolution by a Newton method. Moreover, the use an exact divergence free finite element method velocity-pressure approximation and a discontinuous Galerkin upwind treatment for stresses leads to a robust discretization. A demonstration is provided by the computation of steady solutions at high Weissenberg numbers for the difficult benchmark of the lid driven cavity flow. Numerical results are in good agreement, qualitatively with experiment measurements on real viscoelastic flows, and quantitatively with computations performed by others authors. The numerical algorithm is thus robust. It is also very efficient, as it requires few and mesh-invariant number of linear system resolutions to reach solutions at high Weissenberg number. An adaptive mesh procedure is also proposed, in order to catch accurately both boundary layers and main and secondary vorticies.

We provide new data for future cross validation purpose and point out the existence of a singular point near W e = 2, where the determinant of the Jacobian vanishes and exhibits a change of sign. This singular point has been found quite robust thought mesh refinement. As suggested by both experimental measurements and time-dependent simulations, it could be a bifurcation points associated to a loss of stationarity of the solution. The analysis of this situation by a steady approach is also possible with the tools of the bifurcation theory. In 2009, a pioneering work was performed in this direction by Howell [START_REF] Howell | Dual-mixed finite element approximation of Stokes and nonlinear Stokes problems using trace-free velocity gradients[END_REF] and it should be carried on in the context of the log-conformation formulation.

There are many geometries that could be explored with this approach: contractions, as in the previous reference, and flow around obstacles are interesting benchmarks, while experimental measurements are also available. As new mathematical tools are now available for solving viscoelastic flows problems, the situation becomes mature for considering efficient tridimensional flow simulations. The present log-conformation formulation extends naturally to more complex viscoelastic fluid models, such as Phan-Thien and Tanner, Giesekus or elastoviscoplastic one [START_REF] Saramito | A new constitutive equation for elastoviscoplastic fluid flows[END_REF][START_REF] Saramito | A new elastoviscoplastic model based on the Herschel-Bulkley viscoplasticity[END_REF]. The integration of the Gordon-Schowalter derivative parameter in the present work also open the door of the numerical modeling of liquid crystals, fiber suspensions or active fluids.

W e = 1: 0.9 0.95 1 0.9 0.95 1 software: these expressions contained some errors and, for completeness, the correct expression are provided here. where α 1 = (a -c) 2 + 4b 2 , α 2 = exp((a + c)/2), α 3 = sinh(α 1 /2), and α 4 = cosh(α 1 /2). Notice that the expression in [START_REF] Kane | A comparison of four implementations of the logconformation formulation for viscoelastic fluid flows[END_REF] for exp ′ (χ) should also be fixed in a similar way.

W e = 2:
Observe that the obtained expression degenerates when both a = c and b = 0, i.e. when the matrix is proportional to the identity. In that case χ = a I and exp(χ) = exp(a) I. This possible degeneracy is decided when both |a -c| < ε m and |b| < ε m where ε m stands for the machine epsilon, i.e. the higher positive floating number such that 1 + ε m = 1. Finally, the expressions produced by maxima are simplified and optimized in an automatic way: they are then formatted in fortran, C or C++ languages for a direct use by the Euler-Newton continuation algorithm. The formal calculus approach applies also for a 3 × 3 symmetric matrix: such expression is suitable for tridimensional flow geometries. The explicit expression of the 3×3 matrix exponential is not presented in this paper, as it expands on many pages. Nevertheless, the formal calculus approach is able to automatically generate the corresponding fortran or C++ code in the tridimensional case.

A.2 The κ function and its derivatives

Magnus and Neudecker [43, p. 177] showed that the eigenvalues and eigenvectors are differentiable provided that eigenvalues are simple. Otherwise, when there are multiple eigenvalues, there is a serious problem. Let us study this difficulty by considering the 2 × 2 matrix function:

β(ε, δ) = α + ε δ δ α -ε
where α is some fixed constant and ε and δ are two variables at the vicinity of zero. The two eigenvalues and the two associated unnormalized eigenvectors are

β ± = α ± ε 2 + δ 2 and v ± =   1 - ε δ ± 1 + ε δ 2  
Both eigenvalues are continuous in ε and δ but not differentiable: the conical surface formed by the eigenvalues has a singularity at (ε, δ) = (0, 0) (see Fig. 10.a). For instance

∂β ± ∂ε = ± 1 + δ ε 2 -1 2
For a fixed ratio c = δ/ε however, we can pass through (0, 0) without noticing the singularity, but the derivative depends upon c. Thus, the derivative ∂β ± ∂ε are multi-valued in (0, 0) and the two eigenvalues are not derivable at the origin. Notice that the two eigenvectors only depends upon c.

Nevertheless, the situation is not hopeless for the κ function: Fig. 10.b plots the κ(β, γ) 1,1 component as a function of (ε, δ) for a specific γ value. Other components or γ values show a similar behavior. This observation suggests that κ is differentiable at (ε, δ) = (0, 0) and that its derivative is zero. Let us prove this conjecture. But observing closely the definition (6a)-(6c) of the κ function shows that κ depends only upon (β + -β -)/2 = √ ε 2 + δ 2 and is independent of α. More precisely, from (6b) we have κ1,2 = κ √ ε 2 + δ 2 γ1,2 where with κ(x) = 1 -x/ tanh(x). Then with g(x) = κ′ (x)/x. Observe that lim x→0 g(x) = -2/3. Let us fix the direction c = δ/ε and do ε → 0. By this way v ± depends only on c and so does the unitary matrix q whose columns are the normalized eigenvectors v ± /|v ± |. Thus γ = q T γq does not dependent on ε. Then, the first term of the right-hand side of ( 17) behaves has O(ε). We now turn to the second term of the right-hand side of ( 17): On one hand, we have κ(x) = -x 2 /3 + O(x 3 ) and then κ √ ε 2 + δ 2 = O(ε 2 ). On the other hand, ∂v ± ∂ε = O(δ -1 ) = O(ε -1 ) and so are ∂q ∂ε and ∂γ 1,2 ∂ε . Then, the second term of the right-hand side of ( 17) also behaves has O(ε). Finally, κ1,2 is differentiable in (0, 0) and its derivative is zero. Now, let us turn to κ = qκq T . We have: ∂κ ∂ε = ∂q ∂ε κq T + q ∂κ ∂ε q T + qκ ∂q T ∂ε

∂κ 1,2 ∂ε = ε g ε 2 + δ 2 γ1,2 + κ ε 2 + δ 2 ∂γ 1,2 ∂ε (17)
Recall that q depends only upon the constant direction c. Each term behaves as O(ε) and finally ∂κ ∂ε = 0 at (0, 0). As ε and δ are interchangeable in the expression of κ, a similar deduction leads to ∂κ ∂δ = 0. Remark that, since the derivative is zero, it does not depend on the direction c. To conclude, κ is fully differentiable with respect to β, even when β admits multiple eigenvalues.

The maxima code for computing κ(β, γ) writes: ); q : matrix([v1 [START_REF] Afonso | The finite-volume method in computational rheology[END_REF],v2 [START_REF] Afonso | The finite-volume method in computational rheology[END_REF]], [v1 [START_REF] Amestoy | Multifrontal massively parallel solver (MUMPS). User's guide[END_REF],v2 [START_REF] Amestoy | Multifrontal massively parallel solver (MUMPS). User's guide[END_REF]]); tilde_gamma : transpose(q).gamma.q; hat_kappa(x) := 1 -x/tanh(x); k12 : hat_kappa((beta1-beta2)/2)*tilde_gamma [START_REF] Afonso | The finite-volume method in computational rheology[END_REF][2]; tilde_kappa : matrix([0, k12],[k12,0]); kappa : q.tilde_kappa.transpose(q);

The derivatives of this expression with respect to the matrix coefficients are then treated without any difficulty. As shown in this section, there are two different degenerative cases, when a = c or b = 0. These special cases are decided up to the machine precision and are treated separately. Finally, the expressions produced by maxima are simplified, optimized and formatted in fortran, C or C++ for a direct use by the Euler-Newton continuation algorithm.

Figure 1 :

 1 Figure 1: Incompressible element for the three-field Stokes problem.

Fig. 1

 1 Fig. 1):

Figure 2 :

 2 Figure 2: Anisotropic adaptive mesh.

Figure 3 :

 3 Figure 3: Convergence of the damped Newton algorithm, starting from the solution at W e = 0 for the Oldroyd-B problem on the driven cavity with η p = η v : (a) versus W e for an uniform mesh with h = 1/10 ; (b) versus mesh refinement h for W e = 0.1.

ℵ L 2 h

 2 = 1/20 h = 1/40 h = 1/80 h = 1/160 adapted h min = 10 -2 adapted h min = 5×10 -3

Figure 4 :

 4 Figure 4: The Oldroyd-B problem on the driven cavity with η p = η v . Norm of the solution versus W e.

Fig. 4

 4 Fig.4plots the L 2 norm of the velocity and the log-conformation tensor: observe the good convergence of these quantities with mesh refinement. There are five uniform meshes from h = 1/20 to h = 1/160 and two adapted meshes with h min = 10 -2 and 5 × 10 -3 .

Figure 5 :Figure 6 : 1 ϕ

 561 Figure 5: The Oldroyd-B problem on the driven cavity with η p = η v . Adapted mesh and stream function isovalues for W e = 1 and 2.

Figure 7 : 1 ϕ

 71 Figure 7: The Oldroyd-B problem on the driven cavity with η p = η v . Zoom around vortex for W e = 1.

Figure 8 :Figure 9 :

 89 Figure 8: The Oldroyd-B problem on the driven cavity with η p = η v . Zoom around vortex for W e = 2 and 3.

Figure 10 :

 10 Figure 10: (a) The eigenvalue function β ± = α ± √ ε 2 + δ 2 ; (b) The κ(β, γ) 1,1 function, with γ = ((1, 2), (2, 3)).

  beta : matrix([a,b],[b,c]); gamma : matrix([g00,g01],[g10,g11]); eig : eigenvectors(beta); beta1: eig[1][1][1]; beta2 : eig[1][1][2]; do_unitary(v) := v/sqrt(v[1]**2 + v[2]**2); v1 : do_unitary(eig[2][1][1]); v2 : do_unitary(eig[2][2][START_REF] Afonso | The finite-volume method in computational rheology[END_REF]

Table 1 :

 1 The Oldroyd-B problem on the driven cavity with η p = η v . Maximal and critical Weissenberg number versus mesh refinement: (left) uniform mesh ; (right) adapted mesh.

		3.74	h min W e max W e c
	1/40	2.19	10 -2	3.13	2.33
	1/80 1/160	2.21 2.08	5 × 10 -3	2.78	2.03

Table 2 :

 2 The Oldroyd-B problem on the driven cavity with η p = η v . Comparison of the dimensionless main vortex center position for W e = 1 with others authors. Furthermore, our results are compared quantitatively to previous studies: the locations of the primary vortex center for W e = 1 are listed in table 2. Our results are compared with the results of time-dependent flow using both a finite element method[START_REF] Pan | On the simulation of a time-dependent cavity flow of an Oldroyd-B fluid[END_REF] and a lattice Boltzmann method[START_REF] Su | Lattice Boltzmann method for the simulation of viscoelastic fluid flows over a large range of Weissenberg numbers[END_REF]. It is found that the results of our steady computations are consistent with those of these authors. Figs.7 and.8 show zooms on the left and right secondary vortices. While the main vortex moves from left to right when W e increases, the left vortex grows and the right one decreases in activity. Note that this is also in opposition with inertial effects for viscous Newtonian fluids.

			m	y m			reference
			0.439 0.816		Pan, Hao & Glowinski [51]
			0.433 0.803 Su, Ouyang, Wang, Yang & Zhou [67]
			0.429 0.818			present
	W e	x m	main y m	ϕ min	x l	left y l	ϕ max	x r	right y r	ϕ max
	1 2 3	0.429 0.818 -0.0619 0.0364 0.0388 1.31×10 -6 0.9637 0.0355 9.97×10 -7 0.386 0.828 -0.0555 0.0394 0.0411 1.58×10 -6 0.9663 0.0363 7.07×10 -7 0.335 0.824 -0.0531 0.0485 0.0477 3.22×10 -6 0.9631 0.0347 7.66×10 -7

A Appendices

This appendix contains explicit computations of some complex expressions, namely the exponential of a matrix, the κ function and their derivatives with respect to the matrix coefficients.

A.1 The matrix exponential and its derivatives

The present work uses the free software formal calculus system maxima [START_REF]Maxima, a computer algebra system[END_REF]. With this system, the computation of the exponential of a 2 × 2 symmetric matrix writes: The derivative exp ′ (χ) of this expression with respect to the coefficients of the matrix χ are then treated without any difficulty by this formal calculus system. Notice that Kane, Guénette, and Fortin [35, p. 50] proposed some explicit expressions for exp(χ) and exp ′ (χ), based on the mapple


Cécile Amblard

Stéphane Girard

Ludovic Menneteau

Algebraic properties of copulas defined from matrices

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Algebraic properties of copulas defined from matrices

Cécile Amblard ⋆ , Stéphane Girard, Ludovic Mennetau

We propose a new family of copulas, defined by:

S φ (u, v) = t φ(u)Aφ(v), (u, v) ∈ [0, 1] 2 ,
where φ is a function from [0, 1] to R p and A is a p × p matrix. Let us remark that if p = 2 and A is a diagonal matrix, then S φ reduces to the family proposed in [1]. As a consequence, S φ can be seen as an extension of this former family to arbitrary matrices. First, we shall give sufficient conditions on A and φ to obtain copulas. Then, we shall establish the dependence and symmetry properties of this family of copulas. Finally, we shall study the stability properties of S φ with respect to the operator * (presented for instance in [START_REF] Nelsen | An introduction to copulas[END_REF], p. 194) as well as other algebraic properties.


Yves-Gaël Billet
email: yves-gael.billet@telecom-st-etienne.fr

Christophe Gravier
email: christophe.gravier@telecom-st-etienne.fr

Jacques Fayolle
email: jacques.fayolle@telecom-st-etienne.fr

SWRL-based context awareness for application servers hosting digital services

Keywords: Rules, SWRL, Semantic, context-awareness computing, middleware, digital services, application server

As the number of context-aware applications increases in the real world, it can be quite difficult to deploy such applications in traditional application servers, which are context-agnostics systems. To address this challenge, we propose a novel approach for easing the deployment of contextaware applications into application serversContext is encoded within an OWLdriven knowledge base. We couple this knowledge base with SWRL rules to encode context-awareness thresholds. SWRL rules are not predefined in the application server. They are instead embedded inside the application bundle built by the developer, next to the business logic of the application. At the application deployment time, SWRL rules are extracted to the knowledge base in order to monitor the relevant context for the application to be deployed. At runtime, the context of each session of the application is monitored in the knowledge base. When a rule is triggered (a context-awareness threshold is reached), a broker inside the application server notifies the application so that it adapts its behavior by switching to a more relevant modality. We show how our approach eases the work of developers for building context-aware application by using our context-aware framework.

Needs for a Semantic CAS

We want to provide context awareness for digital service. The idea is to provide a context-aware framework for software architects in order to implement context awareness in their digital services. These services usually run on Applications Servers (henceforth AS) and are consumed through a network using the client-server paradigm. In this situation, context is composed of information about the network and features of the terminal. Context is session specific since applications are hosted on an AS, context-awareness could be externalized as a framework on them.

The objective is to provide a system that supports sensing, perception and reasoning over context. In this way, software developers could only focus on the business logic and adaptation behavior of their application. Behaviors are specific for each application. Software developers must describe and model adaptation behavior for their digital services and communicate it to the context-aware framework.

Semantic Web technologies allow describing context and are comprehensible by both humans and computers and offer a loosely-coupling with the programming language.

Ontologies are used to model domain by describing concepts of the domain and relationships between those concepts [START_REF] Horridge | A Practical Guide To Building OWL Ontologies Using Protege 4 and CO-ODE Tools Edition 1.2[END_REF]. The Web Ontology Language (OWL) was developed to provide a way to represent knowledge understandable by machines using a semantic formalization, in order to facilite computers to interpret human knowledge.

The Semantic Web Rule Language (SWRL) allows us to encode context rules on Horn clauses and designed for OWL. As presented in [START_REF] Horrocks | SWRL: A semantic web rule language combiningOWL and RuleML[END_REF] rules are of the form of an implication between an antecedent (body) and consequent (head). SWRL are logical expression encoded in Conjunctive Normal Forms. It used to classify individuals according conditions. As stated in [START_REF] O'connor | SQWRL: A query language for OWL[END_REF], SPARQL can be used to interrogate an ontology but it is RDF centric so not efficient when using OWL. O'Connor et al. proposes a Semantic Query-enhanced Web Rule Language (SQWRL) based on SWRL.

We propose a semantic context-aware system using OWL, SWRL and SQWRL to describe context and model the adaptation behavior for applications.

Related works

Context-awareness

Among the existing context-aware definitions [START_REF] Zimmermann | An operational definition of context, Modeling and Using Context[END_REF][START_REF] Ak | Understanding and using context, personal and ubiquitous computing[END_REF][START_REF] Strassner | Modelling Context for Autonomic Networking, 5th IEEE International Workshop on Management of Ubiquitous Communications and Services (MUCS)[END_REF], we base our definition from [START_REF] Ak | Understanding and using context, personal and ubiquitous computing[END_REF]. It defines context as "any information that can be used to characterize the situation of an entity. An entity is a person, place, or object that is considered relevant to the interaction between a user and an application, including the user and applications themselves."

We apply this definition in the case of digital services: a context-aware multimedia application is an application (entity), which can use the elements of context (any information that can be used to characterize the situation) to adapt his behavior in order to provide the corresponding interaction to the terminal.

2.2

Architectures of context-aware applications [START_REF] Coutaz | Context is key[END_REF] a global architecture for context-aware applications based on levels of abstraction for a general-purpose context-aware system. They define four levels:

J. Coutaz et al. present in
• Sensing layer provides numeric observables.

• Perception layer provides symbolic observables, which are interpretations from the numeric observable (i.e. transform GPS information in location name).

• Situation and context identification layer identifies context and propose adaptation. It is the reasoning core. • Exploitation layer. It is an adapter between the application and the context infrastructure.

Baldauf et al. in [START_REF] Baldauf | A survey on context-aware systems[END_REF] also use a similar layered architecture to compare middleware and frameworks for context-aware systems.

Existing works on context-aware systems

Current approaches in context-aware frameworks make a clear separation between context acquisition, context processing and use. Thus, to create a separation of concerns as proposed by Dey [START_REF] Ak | Understanding and using context, personal and ubiquitous computing[END_REF]. He proposes a conceptual framework for supporting context-aware applications. The implementation of this framework is known as the Context-Toolkit [START_REF] Ak | Understanding and using context, personal and ubiquitous computing[END_REF]. It aims to ease development and evolution of context-aware applications using an object-based approach. It must be seen as an API for contextaware application. Each object from the API, which can be use to construct an application, has a role (i.e.: collect, transform, aggregate and serve context). The API approach creates a tight coupling between context processing and application.

More recent works, such as CoBrA or SOCAM are infrastructure-based approach rather than API. Hong et al. point out advantages of infrastructure for providing context abstraction [START_REF] Hong | An infrastructure approach to context-aware computing[END_REF].

CoBrA project from Chen et al. [START_REF] Chen | An Intelligent Broker Architecture for Pervasive Context-Aware Systems[END_REF] is agent-oriented. A central unit, so called, a context broker, maintains and manages context on behalf of agents. An agent could be an application running on mobile devices, a service provided by a room, a web service, etc). The broker collects information about context and shares it with agents. This design addresses the issue for providing context-awareness to resource-limited computing devices.

The Service-Oriented Context-Aware Middleware (SOCAM) project introduces by Gu et al. [START_REF] Gu | A middleware for building context-aware mobile services[END_REF] exists as middleware. In SOCAM, different services working together acquire, process, reason and deliver context to agents. The main contribution of this work is the context model, which uses ontology through OWL. The use of ontology allows them to describe context in a semantic way that is independent of programming language.

In context-awareness, the usage is a key element. Works presented in this section are related to the physical world. The main use is to provide services fitted to users activity (e.g.: forward calls to voice mails when users is currently sleeping in the bedroom, switch off light if the room is empty, etc).

Rule driven context-awareness

Our motivations to use a middleware are driven by the will to be unobtrusive [START_REF] Zimmermann | An operational definition of context, Modeling and Using Context[END_REF] and to create an abstraction for context-awareness [START_REF] Ak | Understanding and using context, personal and ubiquitous computing[END_REF]. It aims at deploying regular applications and context-aware applications (related to (1)) and evolving the middleware without impacting the already deployed applications (related to (2)).

We use rules for providing context-awareness thresholds in application through what we call a context signature. It contains rules that define behaviors according to context, for each application. In order to create these rules, the software engineer must define a set of situations. A situation propose a running state for a corresponding environment, as for example, to stream a video through a 4 Mbit/s connection to a terminal we must code the video in 480p.

Each situation is made of observations of the computing environment and a corresponding reaction. Observations describe the computing environment for an application; upon relevant characteristic chosen by the software engineer. Indeed, in our reference scenario, bandwidth is a relevant characteristic. An observation is formatted as a triple (key, comparison operator and value). Values are strings or numeric; in case of a numeric value, a default unit is used. Key is the name of the characteristic, basically, the name of the data gathered by the sensors (e.g.: screen height for an observation about the height of a screen). The operator fills the gap between key and value.

Observations can be as simple as the reference scenario or more complex with additional elements. For example, a situation can be based on observations on bandwidth and screen size. But there can be only one reaction. So, a situation is composed by at least one observable and exactly one consequence.

Each reaction must be a running mode implemented. Rules allow choosing the most appropriate running mode for an application. For each relevant using case of a multimedia application, there must be a situation expressed as a rule. These rules are used against the knowledge base. Once the developer provides the context-aware application as a package which contains the binary and the context signature, the middleware unpacks these elements as shown in Fig. 1.

Fig. 1. Deployment of app in CAAS

The binary is deployed like a standard application and the signature is parsed in order to find rules. As stated before, a rule (situation) is composed with an antecedent (observation) and a consequence (reaction). Each one is then injected as knowledge in the KB. Rules are consider as knowledge because they present running modes and associated necessary conditions for using them, of the newly deployed application.

Context-logic signature defines requirements, in term of context information, for changing application's behavior. We define behavior as a business logic modulation. The main functionality does not change, but can be realized under different forms.

Storing context in a domain ontology

A context aware application server (CAAS) hosts multiple digital services and provides context-awareness for them. Each application provides features under multiple running modes. A user is typified as a session, which consumes features offered by an application. Each session delivers its features through its environment called context. Sensors through context providers can grab it. Each one is able to gather only some kind of information (battery, screen size, localization, etc.) but not all data from the computing environment. In other words, a context-aware application relies on context providers to sense the environment, in order to choose the right business logic modality (or running mode) for a session.

These observations are the foundation for our resulting ontology. It uses two main classes (Context and Session) as shown in Fig. 2.

Context stores all necessary context information about the user computing environment. It is divided into subsets that represent a context provider.

Session represents all active connections to the CAAS. Each user, which consumes an application, is a member of the class session. Like the Session set, this class is divided into subclasses for representing applications.

Relation between a session class and a context class is hasContextInformation. This property is antisymetric, its domain is Session and its range is Context.

Fig. 2. General representation of our ontology

Reference scenario

We will use as reference scenario, a video on demand (VOD) application. The digital service has 3 behaviors (i.e.: 480p, 720p or 1080p) for streaming a video to a user. The context information is the available bandwidth. The application changes the resolution of the video (behavior) according to the user's available bandwidth (context) as shown in Table 1. The context-signature provides semantic formalization of these context configurations that is a set of observations. A context signature maps each context configuration to a specific business logic behavior. This is a loosecoupling solution between business logic and context logic.

We run the above-mentioned ontology against our reference scenario. An application server with our middleware has a context provider called "User Terminal Context" that provide information's about the terminal like bandwidth, screen size and CPU load. The AS hosts a context aware digital service called "Adaptable Videoon-Demand", which provide adaptable video according user's bandwidth. The data gathered by the context provider are represented as datatype relations. In this scenario the context provider gather information about bandwidth, screen size and CPU. As stated before, individuals from the Session context are sorted according to rules, which are injected in the KB. We use SWRL to model the context logic dynamic using this feature. In our reference scenario, the Video-On-Demand service must adapt the video coding according to bandwith. From the OWL point of view, individuals from the set AdaptableVoD must be move either in the c480p, c720p or c1080p subset, each one represents a video coding (480p, 720p and 1080p). This behavior uses the following SWRL:

For more readability, the common part of equations is represents as (0) AdaptableVoD(?s) ^hasContext(?s,?c) ^hasBandwidth( ?c,?b) (0) (0) ^ swrlb:greaterThanOrEqual(?b,3000 ) ^ swrlb:lessThan( ?b,6000) → c480p(?s)

(1)

(0) ^ swrlb:greaterThanOrEqual(?b,6000 ) ^ swrlb:lessThan( ?b,9000) → c720p(?s)

(2)

(0) ^ swrlb:greaterThanOrEqual(?b,9000 ) → c1080p(?s)

(3)

In our model, the reasoner and KB, use OWL and SWRL to choose business logic modality according to context. Each time a new connection is setup, the middleware injects information about the session and the context in the KB. In order to switch modality, the application must be aware of the KB's classification. We use SQWRL to interrogate the ontology in order to notify applications about modality to use for each connection. A typical SQWRL query for this is: Modality(?s) -> sqwrl:select(?s). In our reference scenario: c480p( ?s ) → sqwrl:select(?s) (4) c720p( ?s ) → sqwrl:select(?s)

(5) c1080p( ?s ) → sqwrl:select(?s) [START_REF] Hong | An infrastructure approach to context-aware computing[END_REF] The first one provides all sessions that must use a 480p resolution, the second one for sessions use a 720p and the last one for session uses a 1080p.

Conclusion

We have proposed and implemented a novel architecture that makes AS taking into account the context-awareness of the digital service they host. The architecture employs domain ontology in order to monitor the applications' context. Moreover, when a new context-aware application is deployed on the application server, its context-logic (a set of SWRL rules which rely the domain ontology) is extracted from the application bundle. Each SWRL rule encodes a context-aware threshold for the application. When a rule is triggered at runtime, the application server notifies the application, so that the application change its service delivery modality, as the current context favors another service delivery modality different than the current one.

The primary goal is to help developer of context-aware applications to quickly encode context-aware thresholds. It helps them to develop the context logic of the application for it to adapt its behavior when a significant context change is detected. Developers can seamlessly encode those thresholds by writing the SWRL rules corresponding to the conditions under which the application follows each service delivery modality and then ship the SWRL file into their application bundle.

Previously, this task was not a service offered by the application server, unlike logging, database mapping, authentication, etc., but actually embedded in each application business logics as an ad hoc encoded algorithm. Therefore, this approach is also a framework, as developers of context-aware applications no longer have to write source code for handling context changes.

We are currently developing additional algorithms that will enhance the context assertions in the knowledge base. Especially, we will add to the context logic parameters for gathering context (e.g. context sampling frequency) as each application may have different temporal needs regarding context updates. In future work, we are planning to consider cluster of context-aware application servers.

Table 2 .

 2 Classes of the ontology under our reference scenario

	Name	Parent	Description
	Context Provider	Thing	All context providers
	Session	Thing	All sessions using middleware
	UserTerminalCtx	Context Provider	Context provider that gather information
			about bandwidth, screen size and CPU
			load
	AdaptableVoD	Session	All sessions for the Adaptable Video-
			on-Demand digital service
	C480p	AdaptableVoD	Sessions for AdaptableVoD using 480p
			resolution
	C720p	AdaptableVoD	Sessions for AdaptableVoD using 720p
			resolution
	C1080p	AdaptableVoD	Sessions for AdaptableVoD using
			1080p resolution

Table 3 .

 3 Relations in the ontology under our reference scenario

	Name	Domain	Range	Description
	hasContext	Session	Context	Link a session with a
				context provider
	hasBandwidth	UserTerminalCtx	int	Available bandwidth
				gathered by context
				provider
	hasScreenSize	UserTerminalCtx	string	Screen size gathered by
				context provider
	hasCPULoad	UserTerminalCtx	float	CPU Load gathered by
				context provider