Guy Bayada
email: guy.bayada@insa-lyon.fr

Laurent Chupin
email: laurent.chupin@insa-lyon.fr

Sébastien Martin
email: sebastien.martin@insa-lyon.fr

Viscoelastic fluids in a thin domain

Keywords:

come

Introduction. Much literature research has been devoted to non-Newtonian fluids, in a thin film, in both mathematical aspects and applications. It is well known that numerous biological fluids, blood or physiological secretions like tears or synovial fluids, show these non-Newtonian characteristics. In engineering applications people are interested in controling the flows characteristics to suit various requirements such as maintaining the fluid qualities in a wide range of temperatures and stresses. Introduction of additives lead to non-Newtonian behavior of the modern lubricant. Another application domain is linked to polymers, whose non-Newtonian characteristics appear in a wide range of applications such as the molding or injection processes.

It is to be noticed that, in most practical applications, the geometry of the flow to be considered is anisotropic. A well-known case deals with the study of boundary layers for complex flows [START_REF] Denier | On the boundary-layer equations for power-law fluids[END_REF][START_REF] Fan | Stress boundary layers in the viscoelastic flow past a cylinder in a channel: limiting solutions[END_REF][START_REF] Olagunju | Local similarity solutions for boundary layer flow of a FENE-P fluid[END_REF]. Another case, which is the subject of the present paper, is the lubrication problem in which the fluid is contained between two close surfaces in relative motion. These two applications lead to two very different mathematical models, essentially since the order of magnitude of the parameters in the approximation process is different. For example, in the boundary layer study, the Reynolds number is large and boundary conditions are prescribed at an infinite distance from the solid phase. In lubrication theory, the Reynolds number cannot be too large and boundary conditions are precribed on both surfaces which enclose the fluid. As a consequence, pressure becomes the leading unknown. If such anisotropy can induce some numerical problems in 3D computations, especially as the ratio-aspect of the geometry is sufficiently large, it has however the advantage of allowing some simplification in the equations. So if this approximation process could lead to 2D equations, it could be thought that such simplified equations are easier to solve than the original 3D ones. This explains the amount of work devoted to this topic.

Some particular classes of non-Newtonian models have often been considered. This includes the Bingham flow or the quasi-Newtonian fluids (Carreau's law, the power law or Williamson's law, in which various stress-velocity relations are chosen, see [START_REF] Sac-Épée | On a wide class of non linear models for non-newtonian fluids with mixed boundary conditions in thin domains[END_REF]) and also micropolar ones [START_REF] Bayada | About thin film micropolar asymptotic equations[END_REF]. For these kind of problems, it has been possible to give, in a rigorous way, some thin film approximations of the 3D equations using a so-called generalized Reynolds equation for the pressure. These models, however, considered the fluid as viscous and elasticity effects were neglected. The introdution of such viscoelastic behavior is primilarly described by the Deborah number, denoted De which can be viewed as a measure of the elasticity of the fluid and is related to its relaxation time. One of the laws which seems the most able to describe viscoelastic flows is the Olroyd-B model. This model is based on a constitutive equation which is an interpolation between purely viscous and purely elastic behaviors, thus introducing a supplementary parameter r which describes the relative proportion of both behaviors (the solvant to solute ratio). Considering the Oldroyd model [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF], the momentum, continuity and constitutive equations for an incompressible flow of such a non-Newtonian fluid are, respectively,

ρ ∂ U ∂t ⋆ ⋆ + U ⋆ • ∇U ⋆ -η(1 -r)∆U ⋆ + ∇p ⋆ -div σ ⋆ = 0, (1)
div

U ⋆ = 0, (2)
λ ∂ σ ∂t ⋆ ⋆ + U ⋆ • ∇σ ⋆ + g a (∇U ⋆ , σ ⋆ ) + f (σ ⋆ ) σ ⋆ = 2ηrD(U ⋆ ). (3)
In these equations, ρ, η and λ are positive constants which respectively correspond to the fluid density, the fluid viscosity and the relaxation time. Equations ( 1)-(3) make up a system of 10 scalar equations with 10 unknowns: the lubricant velocity vector U ⋆ = (u ⋆ 1 , u ⋆ 2 , w ⋆ ), the pressure p ⋆ and the extra-stress symmeric tensor σ ⋆ = (σ ⋆ i,j ) 1≤i,j≤3 . The bilinear application g a , -1 ≤ a ≤ 1, is defined by

g a (∇U ⋆ , σ ⋆ ) = σ ⋆ • W (U ⋆ ) -W (U ⋆ ) • σ ⋆ -a(σ ⋆ • D(U ⋆ ) + D(U ⋆ ) • σ ⋆ )
where D(U ⋆ ) and W (U ⋆ ) are respectively the symmetric and skew-symmetric parts of the velocity gradient ∇U ⋆ . Usually, D(U ⋆ ) is called the rate of strain tensor and W (U ⋆ ) is called the vorticity tensor. Notice that the parameter a is considered to interpolate between upper convected (a = 1) and lower convective derivatives (a = -1), the case a = 0 being the corotational case [START_REF] Joseph | Fluid dynamics of viscoelastic liquids[END_REF]. Note that taking r = 1 allows us to recover various forms of the generalized Maxwell model. Then, by choosing f as the identity, this model is the classical Maxwell one. By introducing a linearized form of f (see in particular [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF]), the Phan Tien-Tanner laws [START_REF] Tanner | Rheology: an historical perspective[END_REF] are obtained. Conversely, a Newtonian flow is described by choosing r = 0.

From the mathematical point of view, few results exist concerning the existence or uniqueness of a solution for true 3D or 2D viscoelastic models [START_REF] Chupin | Some theoretical results concerning diphasic viscoelastic flows of the Oldroyd kind[END_REF][START_REF] Guillopé | Existence results for the flow of viscoelastic fluids with a differential constitutive law[END_REF][START_REF] Molinet | On the global and periodic regular flows of viscoelastic fluids with a differential constitutive law[END_REF], also the way to obtain the related thin film approximation is mainly heuristic. A primary approach, which is often used in engineering literature, is to take the parameter defining the (relative) thickness of the flow as the leading small parameter and to use the Deborah number as a pertubation parameter. This has been carried out in the lubrication field by Tichy [START_REF] Tichy | Non-Newtonian lubrication with the convective Maxwell model[END_REF] starting from the upper convected Maxwell model (r = 1, f = Id, a = 1). The case of a Deborah number of the same order of magnitude as the relative thickness has been studied by Tichy and Huang for the UCM Maxwell model and by Bellout [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF] for the Phan Tien-Tanner model. In all these researchs, a nonlinear Reynolds equation is gained, allowing the pressure in the thin film to be directly computed. The same procedures can also include the free boundary upper surface of the flow (thin coating problem) or the inertia [START_REF] Khayat | Transient free-surface flow of a viscoelastic fluid in a narrow channel[END_REF][START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF][START_REF] Zhang | Surfactant spreading on a thin weakly viscoelastic film[END_REF]. However, the goal of these last studies is different, as the primary unknown is not an equation for the pressure but an equation describing the evolution of the free boundary (a generalized shallow water equation).

The present paper addresses the mathematical and numerical study of a large class of viscoelastic thin film flows described by an Olroyd-B model in which the Deborah number has the same order of magnitude as the thickness of the fluid. This assumption allows the order of Newtonian and non-Newtonian contributions (see [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF] for mechanical comments) to be balanced. Boundary conditions are chosen to be applied to the usual lubrication problems. After scaling both equations and the stress tensor in an adequate way, we are able to obtain an asymptotic 2D problem. This problem generalizes the work of Bellout and Tichy, and concerns not only the rheological model but also can take the 2D dimension (instead of 1D for the pressure asymptotic problem) into account. Obtaining the asymptotic problem is partly an heuristic process, so we have to rigorously prove the solvability of this problem. This is the goal of Section 2 which is divided in two parts for sake of clarity. The newtonian case (r = 0) is studied first and a new way to obtain an existence and uniqueness result for the problem is proposed using velocity as a leading unknown. This type of approach can be easily generalized to the viscoelastic case by using a monotonicity property of the nonlinear term. Interestingly, an existence and uniqueness result is obtained exactly for the same range of the r parameters as in the initial 3D problem. In numerous problems in thin fields, it is possible to eliminate the velocity in the limit problem, thus only retaining a Reynolds equation with respect to the pressure. It is different in our case and we have to solve a nonlinear coupled problem in which a degenerate Stokes equation is still present. A new algorithm related to the Uzawa one is presented and the convergence theorems are given. Lastly, numerical comparisons between various models are given and the importance of obtained 2D and not only a 1D approximation is emphasized.

1. Mathematical formulation. The space coordinates are denoted by

(x ⋆ 1 , x ⋆ 2 , z ⋆ ) or more simply by (x ⋆ , z ⋆ ) with x ⋆ = (x ⋆ 1 , x ⋆ 2 ).
Let ω be a fixed bounded domain of the plane z ⋆ = 0. We suppose that ω has a Lipschitz continuous boundary ∂ω. The upper surface of the gap is defined by z ⋆ = H(x ⋆ ) with H ∈ C 1 (ω). Let Ω be the following set (see Fig. 1):

Ω = {(x ⋆ , z ⋆ ) ∈ R 3 , x ⋆ ∈ ω and 0 < z ⋆ < H(x ⋆ )}. z = H(x 1 , x 2 ) ∂ω ω z 0 x 2 x 1 Ω Fig. 1.
The physical domain 1.1. Thin film flow equations. By introducing characteristic lengths L for the domain ω and H for the size of the gap, we can define the ratio

ε = H L
which is, in the realistic physical case of lubrification, very small. The governing equations ( 1)-( 3) can be expressed in dimensionless form in terms of the following dimensionless quantities :

x = x ⋆ L , z = z ⋆ εL , u i = u ⋆ i U , w = w ⋆ εU , (4)
p = p ⋆ ε 2 L η U , σ = σ ⋆ εL ηU , t = t ⋆ U L . (5)
We now introduce two classical numbers : the Reynolds number Re which characterises the viscous forces compared to the convective ones and the Deborah number De which highlights the elasticity of the fluid. They are defined by

Re = ρUL η , De = λU εL . (6)
For convenience, we also introduce the normalized gap function:

h(x) = H(x) εL . (7)
Remark 1.1. This scaling process is motivated by the following considerations:

• The length and velocity scaling (4) takes into account the thin film nature of the lubrication flow. • Classically, in lubrication theory, if the horizontal shear velocity is of the order of 1, then the real pressure is of the order of 1/ε 2 . Such result can be obtained by proving that ε 2 p is bounded in a suitable norm (see [START_REF] Bayada | The transition between the Stokes equations and the Reynolds equation: a mathematical proof[END_REF] for a rigorous mathematical explanation). • If we want to balance the Newtonian and non-Newtonian contribution, we must assume that the stress tensor is of the order of 1/ε and the usual Deborah number λU/L is of the order of ε (see [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF] for further explanations).

By substituting these dimensionless variables (4)- [START_REF] Denier | On the boundary-layer equations for power-law fluids[END_REF] in Equations ( 1)-( 3) we obtain the dimensionless governing equations.

• The three components of the momentum equation ( 1) are written

                                             Re d u 1 dt -(1 -r) ∂ 2 u 1 ∂x 2 1 + ∂ 2 u 1 ∂x 2 2 + 1 ε 2 ∂ 2 u 1 ∂z 2 + 1 ε 2 ∂ p ∂x 1 - 1 ε ∂ σ 1,1 ∂x 1 + ∂ σ 1,2 ∂x 2 + 1 ε ∂ σ 1,3 ∂z = 0, Re d u 2 dt -(1 -r) ∂ 2 u 2 ∂x 2 1 + ∂ 2 u 2 ∂x 2 2 + 1 ε 2 ∂ 2 u 2 ∂z 2 + 1 ε 2 ∂ p ∂x 2 - 1 ε ∂ σ 1,2 ∂x 1 + ∂ σ 2,2 ∂x 2 + 1 ε ∂ σ 2,3 ∂z = 0, εRe d w dt -ε(1 -r) ∂ 2 w ∂x 2 1 + ∂ 2 w ∂x 2 2 + 1 ε 2 ∂ 2 w ∂z 2 + 1 ε 3 ∂ p ∂z - 1 ε ∂ σ 1,3 ∂x 1 + ∂ σ 2,3 ∂x 2 + 1 ε ∂ σ 3,3 ∂z = 0.
When ε tends to zero, these equations formally reduce to the following set of equations:

               -(1 -r) ∂ 2 u 1 ∂z 2 + ∂ p ∂x 1 - ∂ σ 1,3 ∂z = 0, -(1 -r) ∂ 2 u 2 ∂z 2 + ∂ p ∂x 2 - ∂ σ 2,3 ∂z = 0, ∂ p ∂z = 0, (8)
• Due to the previous dimensionless procedure the free divergence condition is preserved for the dimensionless variables:

∂ u 1 ∂x 1 + ∂ u 2 ∂x 2 + ∂ w ∂z = 0, (9)
• Concerning the constitutive law, the process is similar: equations are written for the dimensionless quantities, then, passing formally to the limit ε → 0, the following equations are obtained:

                                         σ 1,1 + De(1 -a)σ 1,3 ∂ u 1 ∂z = 0, σ 2,2 + De(1 -a)σ 2,3 ∂ u 2 ∂z = 0, σ 3,3 -De(1 + a)(σ 1,3 ∂ u 1 ∂z + σ 2,3 ∂ u 2 ∂z ) = 0, σ 1,2 + De 2 (1 -a)(σ 2,3 ∂ u 1 ∂z + σ 1,3 ∂ u 2 ∂z ) = 0, σ 1,3 + De 2 (1 -a)σ 3,3 ∂ u 1 ∂z -(1 + a)σ 1,2 ∂ u 2 ∂z -(1 + a)σ 1,1 ∂ u 1 ∂z = r ∂ u 1 ∂z , σ 2,3 + De 2 (1 -a)σ 3,3 ∂ u 2 ∂z -(1 + a)σ 1,2 ∂ u 1 ∂z -(1 + a)σ 2,2 ∂ u 2 ∂z = r ∂ u 2 ∂z . (10)
In this system, it is easy to see that coefficients σ 1,1 , σ 2,2 , σ 3,3 and σ 1,2 can be expressed as a fonction of σ 1,3 , σ 2,3 and of the velocity (u 1 , u 2 ). In addition, using the last two equations, σ 1,3 and σ 2,3 are expressed with respect to the velocity:

σ 1,3 = r ∂ u 1 ∂z 1 + De 2 (1 -a 2 ) ∂ u 1 ∂z 2 + ∂ u 2 ∂z 2 , σ 2,3 = r ∂ u 2 ∂z 1 + De 2 (1 -a 2 ) ∂ u 1 ∂z 2 + ∂ u 2 ∂z 2 .
For the sake of simplicity, let us denote by u the first two coordinates of the velocity vector : u = (u 1 , u 2 ) and by β the following two components of the stress tensor : β = (σ 1,3 , σ 2,3 ). The system obtained can be written in the following form:

                     -(1 -r) ∂ 2 u ∂z 2 - ∂ β ∂z + ∇ x p = 0, with β = r ∂ u ∂z 1 + De 2 (1 -a 2 ) ∂ u ∂z 2 , ∂ p ∂z = 0, div x u + ∂ w ∂z = 0, (11)
all the other components of the stress tensor being directly deduced from equations [START_REF] Khayat | Transient free-surface flow of a viscoelastic fluid in a narrow channel[END_REF].

The vertical velocity w can be deduced from the horizontal velocity u by the free divergence condition. More clearly, when taking the boundary conditions [START_REF] Tanner | Rheology: an historical perspective[END_REF] presented in the next subsection into account, we can check that problem [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] implies

                       -(1 -r) ∂ 2 u ∂z 2 -r ∂ ∂z    ∂ u ∂z 1 + De 2 (1 -a 2 ) ∂ u ∂z 2    + ∇ x p = 0, ∂ p ∂z = 0, div x h 0 u dz = w(•, 0) -w(•, h). (12)
1.2. Boundary conditions. System (12) will be the subject of the forthcoming theoretical study, as it enables the knowledge of the pressure p (the primary factor of interest in lubrication problems) and the horizontal velocity u to be obtained (while the vertical velocity w is in the real variables of the order of ε). Let us now introduce the boundary conditions. As it is well-known (see [START_REF] Assemien | Inertial effects in the asymptotic behavior of a thin film flow[END_REF]), passing from 3D problems to 2D ones may cause boundary layer phenomena on the lateral parts of Ω. So only a part of the boundary condition for the initial problem has to be considered in the study of [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF]. We have to retain the following typical (no-slip) boundary conditions at z = 0 and z = h:

• u(•, 0) = s and u(•, h) = 0 on ω,

• w(•, 0) = 0 and w(•, h) = 0 on ω, in which s := (s 1 , s 2 ) is a given velocity vector. Moreover, two kinds of boundary conditions can be considered along this lateral boundary, one associated to the pressure data, the other one to the average flux data. The choice of the conditions depends highly on the devices to be considered. In most of the physical problems, two types of boundary conditions are simultaneously used: Neumann-type conditions and Dirichlet conditions. Thus, in the general case, the set of equations [START_REF] Lions | Variational inequalities[END_REF] has to be considered with the following boundary conditions:

p = p 0 on ∂ω p , h 0 u dz • n = q 0 on ∂ω q , ( 13
)
where ∂ω p and ∂ω q define a partition of the boundary ∂ω, and n denotes the outward normal vector on ∂ω. Notice that ∂ω p (resp. ∂ω q ) may be the union of a finite number of connected components denoted ∂ω p i (resp. ∂ω q i ) (see Fig. 2). Let us notice that a compatibility condition on the total flux is needed if ∂ω p = ∅: 11) has a non classical feature due to the nonlinear elastic term. We propose in this section a two-step procedure. In the first one, the nonlinear term is cancelled, thus obtaining a purely Newtonian model. The idea is to solve this problem in terms of velocity first. This is completely different from the usual approach [START_REF] Bayada | The transition between the Stokes equations and the Reynolds equation: a mathematical proof[END_REF] in which the velocity is expressed in terms of the pressure to obtain a problem in which the pressure is the only unknown. The advantages of this new procedure clearly appear in the second step as introducing the nonlinear term (to cope with viscoelastic effects) is then straightforward.

∂ω q 0 = 0. ∂ω q 1 ∂ω p 1 ∂ω q 2 ∂ω p 2 ω Pressure imposed

Flux imposed

2.1. The newtonian case. The Newtonian case corresponds to the case where the stress tensor σ is zero. In the limit equations [START_REF] Lions | Variational inequalities[END_REF], this means that β = 0. In this subsection, we first state the strong and weak formulations of the problem. So we do not only provide a rigorous mathematical study but also establish the relevance of the weak formulation with respect to the physical (strong) formulation. Thus, let us introduce the formulations in the purely Newtonian case.

Strong formulation:

The problem deals with boundary conditions of two types: Neumann conditions and non-homogeneous Dirichlet conditions. By introducing some kind of source-term, it is possible to obtain an equivalent problem with homogeneous Dirichlet conditions. In fact let p 0 be an extension of p 0 on the closed set ω. It is obviously equivalent of working with a reduced pressure p = pp 0 instead of the effective pressure p. Due to this change the strong formulation (P s ) is slightly modified by the introduction of a non-zero right-hand side F = -∇ x p 0 (instead of 0), which takes into account the translation of the pressure. In the whole study, we will consider that the following assumptions on the data hold:

Assumption 1 (Regularity of the data). • h ∈ C 0 (ω), h ≥ h 0 > 0, • s ∈ L 2 (Γ -), where Γ -(resp. Γ + ) denotes the lower (resp. upper) boundary of Ω, i.e. Γ -= {(x, 0), x ∈ ω}, Γ + = {(x, h(x)), x ∈ ω}. • F ∈ L 2 (Ω), • q 0 ∈ L 2 (∂ω q ).
Now, the strong formulation is as follows:

(P s )                                            - ∂ 2 u ∂z 2 + ∇ x p = F , in L 2 (Ω), (14)
∂ p ∂z = 0, in L 2 (Ω), (15)
div x h 0 u(•, z) dz = 0, in L 2 (ω), (16)
u = s, in L 2 (Γ -), (17)
u = 0, in L 2 (Γ + ), (18)
p = 0, in L 2 (∂ω p ), (19)
h 0 u(•, z) dz • n = q 0 , in L 2 (∂ω q ). ( 20
)
It can be noticed that this set of equations can be reduced to the classical Reynolds equation (see in particular [START_REF] Bayada | The transition between the Stokes equations and the Reynolds equation: a mathematical proof[END_REF]). In fact, integrating Equation ( 14) twice with respect to z (and taking into account the velocity boundary conditions ( 17)-( 18)), the velocity u is obtained as a function of the pressure p. Then, putting this expression into Equation ( 16) gives:

div h 3 6 ∇p = div (sh) .
In this purely Newtonian case, the Reynolds formulation enables a straightforward existence and uniqueness result (via elliptic theory) to be given. In this paper we propose an alternate approach which will be easily adapted to the viscoelastic case (although the Reynolds approach could not be easily extended to this nonlinear case).

Weak formulation:

First, let us introduce the functional space which is used in the weak formulation. For s ∈ R 2 and q 0 ∈ L 1 (∂ω q ), we define the following subspace of L 2 (Ω):

K(s, q 0 ) = ϕ ∈ L 2 (Ω), ∂ ϕ ∂z ∈ L 2 (Ω), ϕ = s in L 2 (Γ -), ϕ = 0 in L 2 (Γ + ), ∀θ ∈ {θ ∈ D(ω) s.t. ∃ζ ∈ R, θ| ∂ω p = ζ}, ω ∇ x θ • h 0 ϕ(•, z) dz = ∂ω q (θ -ζ) q 0 .
The space K(s, q 0 ) is equipped with the norm:

ϕ z = Ω ∂ ϕ ∂z 2 1/2 . ( 21
)
Remark 2.1. It can be noticed that for every function ϕ satisfying

ϕ ∈ L 2 (Ω), ∂ ϕ ∂z ∈ L 2 (Ω),
it is possible to define, by density, its trace γ

-(ϕ) ∈ L 2 (Γ -) on Γ -(resp. γ + (ϕ) ∈ L 2 (Γ + )).
Thus, the boundary conditions on Γ -and Γ + in the definition of K(s, q 0 ) make sense. However, the indexes γ ± will be dropped for the sake of simplicity.

The following proposition will enable us to better understand the interest of this set:

Proposition 2.1. Let ϕ ∈ H 1 (Ω). We have the following equivalence ϕ ∈ K(s, q 0 ) ⇐⇒ ϕ satifies Equations ( 16)-( 18) and [START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF].

Proof. It is clear that if ϕ ∈ K(s, q 0 ), then Equations ( 17)-( 18) hold (see the definition of the functional space). Now, using an integration by parts, if ϕ ∈ K(s, q 0 ) ∩ H 1 (Ω) then for all θ ∈ D(ω), θ being constant on ∂ω p , we obtain

- ω θ div x h 0 ϕ(•, z) dz + ∂ω q θ h 0 ϕ(•, z) dz • n = ∂ω q θ q 0 .
In particular, for all θ ∈ D(ω), we find

ω θ div x h 0 ϕ(•, z) dz = 0
thus Equation ( 16) holds. Then, for all θ ∈ D(∂ω q ), extended on ω such that θ ∈ D(ω) and θ| ∂ω p = 0, we obtain [START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF] holds. This concludes the proof of the necessary condition. This condition is clearly sufficient.

∂ω q θ h 0 ϕ(•, z) dz • n -q 0 = 0, i.e. Equation
The method of attaining the weak formulation of the problem will be now provided. Let (u, p) be a regular solution of ( 14)- [START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF], and let ϕ ∈ K(s, q 0 ). Multiplying Equation ( 14) by uϕ and integrating over Ω, we obtain

Ω - ∂ 2 u ∂z 2 • (u -ϕ) + Ω ∇ x p • (u -ϕ) = Ω F • (u -ϕ).
Since uϕ ∈ K(0, 0), we can integrate by parts the first integral, and use p as a test function like θ to cancel the second integral (let us recall here that p does not depend on z). In particular, we deduce the weak formulation of the problem:

(P w )    Find u ∈ K(s, q 0 ) such that Ω ∂ u ∂z • ∂ (u -ϕ) ∂z ≤ Ω F • (u -ϕ), ∀ ϕ ∈ K(s, q 0 ). ( 22
)
Now, this subsection is concluded with two major results. We first give an existence and uniqueness result for the weak problem and then we describe the link between the two formulations.

Theorem 1 (Newtonian case). Problem (P w ) admits a unique solution.

Proof. The proof is based on the theory of variational inequalities [START_REF] Lions | Variational inequalities[END_REF]. Obviously, the space (K(s, q 0 ), • z ) is closed in L 2 (Ω). Moreover linearity of the boundary conditions leads to the affine property of the space so that it is convex. Thus it remains to prove that the space is non-empty. Using Proposition 2.1, we look for a function satisfying Equations ( 16)-( 18) and [START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF]. It is obvious that the function

φ = a 2 z(z -h)) + s h -z h , (23)
satisfies [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF] and [START_REF] Tanner | Rheology: an historical perspective[END_REF]. Here, a is any vector which only depends on x (to be further detailed). In order to ensure that φ satisfies Equations ( 17), ( 18) and ( 20), a has to satisfy:

       div h 3 12 a = div sh 2 on ω, sh 2 - h 3 12 a • n = q 0 on ∂ω q . (24)
In order to state that a certain a, satisfying the earlier set of equations exists, we consider the following Reynolds problem (as an auxiliary problem):

       div h 3 12 ∇π = div sh 2 on ω, sh 2 - h 3 12 ∇π • n = q 0 on ∂ω. ( 25
)
Obviously a unique π ∈ H 1 (ω)/R satisfying (25) exists. So, choosing a = ∇π ∈ L 2 (ω), the proof is concluded as by definition, a satisfies Equations ( 17), [START_REF] Tanner | Rheology: an historical perspective[END_REF] and [START_REF] Zhang | Non-Newtonian effects on lubricant thin film flows[END_REF]. Thus, the function φ defined by Equation (23), with the previous choice for a, belongs to K(s, q 0 ) which is consequently non-empty.

The link between (P w ) and (P s ) is given by the following theorem.

Theorem 2 (Newtonian case). Let u be the unique solution of (P w ). (i) A unique p ∈ H 1 (ω) such that (u, p) satisfies ( 14), ( 15), ( 17)-( 19) exists.

(ii) Moreover, if u ∈ H 1 (Ω), then ( 16) and ( 20) hold. In particular, (u, p) is the unique solution of (P s ).

Proof. The result is checked in three steps:

• Step 1: Let us state that Equations ( 14) and ( 15) hold.

For this, we use the de Rham theorem in order to ensure the existence of a pressure p. Choosing ϕ = u ± ϕ with ϕ ∈ K(0, 0) ∩ D(Ω) as a test function in Equation ( 22), we deduce that

∀ϕ ∈ K(0, 0) ∩ D(Ω), Ω ∂ u ∂z • ∂ ϕ ∂z = Ω F • ϕ. (26)
Then, as u belongs to K(s, q 0 ), we find

∀ϕ ∈ K(0, 0) ∩ D(Ω), - ∂ 2 u ∂z 2 -F , ϕ = 0, ( 27
)
where •, • denotes the classical duality product. The next lemma allows us to use the classical De Rham theorem to find a pressure:

Lemma 2.1. For ϕ = (ϕ 1 , ϕ 2 ) ∈ K(0, 0)∩D(Ω), ϕ 3 ∈ D(Ω) exists such that div(ϕ 1 , ϕ 2 , ϕ 3 ) = 0. Conversely, if Φ = (ϕ 1 , ϕ 2 , ϕ 3 ) ∈ D(Ω) is such that div Φ = 0, then (ϕ 1 , ϕ 2 ) ∈ K(0, 0). Proof. For ϕ ∈ K(0, 0) ∩ D(Ω) it is sufficient to define ϕ 3 (x, z) = - z 0 div x ϕ(x, ξ) dξ so that div x ϕ + ∂ ϕ 3 ∂z = 0, ( 28
)
with ϕ 3 ∈ D(Ω). Conversely let us write ϕ = (ϕ 1 , ϕ 2 ). If Equation (28) holds, then using the fact that ϕ 3 is zero at the boundaries z = 0 and z = h, we find

div x h 0 ϕ(•, z) dz = 0. Moreover, if ϕ ∈ D(Ω) then ϕ = 0, on ∂ω q h 0 ϕ(•, z) dz • n = 0, on ∂ω q , i.e. ϕ ∈ K(0, 0). Let us define F = - ∂ 2 u ∂z 2 + ∇ x p 0 , 0 ,
and using the previous lemma, Equation ( 27) is rewritten as:

∀Φ ∈ D(Ω) such that div Φ = 0, F , Φ = 0.
With the De Rham theorem, we deduce that a unique pressure p ∈ D ′ (Ω)/R exists such that F = ∇p, with

∇p = ∇ x p, ∂ p ∂z , so that - ∂ 2 u ∂z 2 + ∇ x p = F , in D ′ (Ω), (29)
∂ p ∂z = 0, in D ′ (Ω), (30)
Now let us discuss the regularity of u and p: as u is a solution of problem (P w ), then

u, ∂ u ∂z ∈ L 2 (Ω).
In particular, if u is extended by

0 on {(x, z) ∈ Ω, z ≥ h(x)}, by denoting h ∞ = h L ∞ (ω) , we obtain u ∈ C([0, h ∞ ]; L 2 (ω)).
Now, by Equation (29), as F and p do not depend on z, we obtain

∂ ∂z ∂ 2 u ∂z 2 ∈ L 2 ((0, h ∞ ); L 2 (ω)),
and since

∂ 2 u ∂z 2 ∈ D ′ ((0, h ∞ ); L 2 (ω)), then ∂ 2 u ∂z 2 ∈ C([0, h ∞ ]; L 2 (ω)).
Moreover, by Equation (29),

∇p ∈ C([0, h ∞ ]; L 2 (ω))
and by Equation (30) (p does not depend on z), we conclude that ∇p ∈ L 2 (ω), i.e. (u, p) satisfies ( 14) and [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF]. In particular, the boundary conditions for the pressure on ∂ω make sense.

• Step 2: Let us state that Equation [START_REF] Tichy | Non-Newtonian lubrication with the convective Maxwell model[END_REF] holds.

The only point to be checked consists in showing that the pressure p is constant along the curve ∂ω p . Since u is a weak solution (that is to say a solution of (27)) and (u, p) satisfies ( 14), we immediately deduce by difference that the pressure p satisfies:

∀ϕ ∈ K(0, 0) ∩ D(Ω) Ω ∇ x p • ϕ = 0. ( 31
)
The end of this section thus will be devoted to show that this condition (31) implies that p is constant on ∂ω p . The proof is realised in three sub-steps:

⊲ Step 2-1. (Technical lemma)
Lemma 2.2. The following application is surjective:

Φ : K(0, 0) ∩ D(Ω) -→ X = {f ∈ D(ω), div x f = 0, f • n = 0 on ∂ω q } ϕ -→ h 0 ϕ(•, z) dz
Proof. Using Proposition 2.1, we show that this application is well defined and with values in X. For f ∈ X, we define Then, for ϕ ∈ D(∂ω p ), extended by zero on ∂ω q , we have ϕ + ϕ ∈ Z, so that

ϕ(x, z) = 12 h(x) 3 z(z -h(x))f (x)
∂ω p ∂(ϕ + ϕ) ∂τ = 0.
By making the difference with Equation (32), we find

∀ ϕ ∈ D(∂ω p ), ∂ω p p ∂ ϕ ∂τ = 0,
that is to say p is constant along each connected component of ∂ω p .

⊲

Step 2-3. (The value of the boundary pressure is the same on each connected component of ∂ω p ) -If ∂ω p is composed of n connected components ∂ω p i , i ∈ {1, ..., n}, then its complementary subset is also composed of n connected components: ∂ω q i , i ∈ {1, ..., n} (see Fig. 2). For all (a 1 , ..., a n ) ∈ R n , we define a function a ∈ D(∂ω) such that a| ∂ω q i = a i . We obtain ϕ + ϕ + a ∈ Z, so that

∂ω p ∂(ϕ + ϕ + a) ∂τ = 0.
By making the difference between (2.1) and (2.1) we find ∂ω p ∂a ∂τ = 0.

Since p is constant on each connected component ∂ω p i of ∂ω p (with the value p i ), this equality may be also written as

0 = ∂ω p ∂a ∂τ = ∂ω p p ∂a ∂τ = n i=1 p i ∂ω p i ∂a ∂τ = n i=1 p i (a i -a i-1 ) = n i=1 a i (p i -p i+1 ), (33)
with the convention p 1 = p n+1 and a 0 = a n . As Equation (33) must be satisfied for all constants a i ∈ R, we find that all the p i have the same value.

• Step 3: Since u ∈ K(s, q 0 ), Equations ( 17) and ( 18) hold. If furthermore u ∈ H 1 (Ω) then, by Proposition 2.1, Equations ( 16) and ( 20) hold.

2.2. The viscoelastic case. The introduction of viscoelastic phenomena differs from the purely Newtonian case by the effect of nonlinear additive terms. However, we show in this subsection that the approach developed earlier allows us to analyse rigorously the complete problem. Due to the introduction of the nonlinear terms, the mathematical analysis of the corresponding weak formulation (to be detailed further) has to be adapted in order to ensure the existence and uniqueness of the (weak) solution. Thus, let us first introduce the strong and weak formulations of the viscoelastic problem in a thin domain.

Strong formulation:

(Q s )                                                    -(1 -r) ∂ 2 u ∂z 2 -r ∂ ∂z    ∂ u ∂z 1 + C 2 ∂ u ∂z 2    + ∇ x p = F , in L 2 (Ω), (34)
∂ p ∂z = 0, in L 2 (Ω), (35)
div x h 0 u(•, z) dz = 0, in L 2 (ω), ( 36
) u = s, in L 2 (Γ -), ( 37
) u = 0, in L 2 (Γ + ), ( 38
) p = 0, in L 2 (∂ω p ), ( 39
) h 0 u(•, z) dz • n = q 0 , in L 2 (∂ω q ). ( 40
)
where the constant C ≥ 0 includes viscoelastic parameters, namely C 2 = De 2 (1a 2 ).

Weak formulation:

Nonlinear terms due to viscoelasticity have to be taken into account, leading to a significant modification of the Newtonian case, so that the weak formulation of the problem is written as:

(Q w ) Find u ∈ K(s, q 0 ) such that ≪ Au, u -ϕ ≫ ≤ ≪ F , u -ϕ ≫, ∀ ϕ ∈ K(s, q 0 ), (41)
where A : K(s, q 0 ) → (K(s, q 0 )) ′ is the operator defined by

≪ Au, v ≫= (1 -r) ∂ u ∂z , ∂ v ∂z L 2 (Ω) + r    ∂ u ∂z 1 + C 2 ∂ u ∂z 2 , ∂ v ∂z    L 2 (Ω)
and let us recall that K(s, q 0 ) is equipped with the norm • z (see its definition given by ( 21)). Now we give the following theorem, which is a generalisation of Theorem 1 taking into account the viscoelastic terms.

Theorem 3 (Viscoelastic case). If r < 8/9, problem (Q w ) admits a unique solution.

Proof. The proof is based on a classical result on variational inequalities with monotone operators (see [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF], page 247). It is obtained using three steps:

• Step 1: bounded operator. Obviously, since r ≥ 0, we write

≪ Au, u ≫= (1 -r) Ω ∂ u ∂z 2 + r Ω ∂ u ∂z 2 1 + C 2 ∂ u ∂z 2 ≤ u 2 z ,
which means that A is bounded.

• Step 2: coercive operator. Here, we use the fact that r < 1: indeed,

≪ Au, u ≫ ||u|| z ≥ (1 -r)||u|| z , so that lim u z →+∞ ≪ Au, u ≫ ||u|| z = +∞.
• Step 3: monotone operator. We show here that the operator A is strictly monotone if and only if r < 8/9 (independently of the constant C). Thus let us focus on the formula

≪ Au -Av, u -v ≫: ≪ Au -Av, u -v ≫ = (1 -r) Ω ∂ u ∂z - ∂ v ∂z 2 + r Ω ∂ u ∂z - ∂ v ∂z 2 1 -C 2 ∂ u ∂z • ∂ v ∂z 1 + C 2 ∂ u ∂z 2 1 + C 2 ∂ v ∂z 2 = Ω ∂ u ∂z - ∂ v ∂z 2 B ∂ u ∂z , ∂ v ∂z 1 + C 2 ∂ u ∂z 2 1 + C 2 ∂ v ∂z 2 where B(a, b) = (1 -r) 1 + C 2 a 2 1 + C 2 b 2 + r 1 -C 2 a • b . Rewriting this term as B(a, b) = C 2 (1 -r) 1 + C 2 b 2 a - rb 2(1 -r) 1 + C 2 b 2 2 + 1 4(1 -r) 1 + C 2 b 2 2(1 -r)C 2 b 2 + 3r 2 -12r + 8 4(1 -r) 2 + r 3 64(1 -r) 3 1 + C 2 b 2 (8 -9r),
we deduce the sign of ≪ Au -Av, uv ≫. In fact, studying the sign of B(a, b) gives the following:

⊲ if r < 8/9
, the operator A is strictly monotone. ⊲ if r = 8/9, the operator A is monotone. ⊲ if r > 8/9, the operator A is non monotone: we can find u and v such that

≪ Au -Av, u -v ≫< 0.
Now, the proof is concluded using the theory of monotone operators in variational inequalities.

Remark 2.2 (A non-uniqueness result). Interestingly, we can prove that the problem is well-posed if r < 8/9. When r = 8/9, the proof of Theorem 3 ensures the existence of a weak solution (but not necessarily uniqueness). When r > 8/9, thus does not even state an existence result. However, using a simple geometrical configuration (h ≡ 1), a counter-example for uniqueness can be established (see [START_REF] Boyer | Numerical study of viscoelastic mixtures through a Cahn-Hilliard flow model[END_REF] for further details).

The link between (Q w ) and (Q s ) is given by the following theorem.

Theorem 4 (Viscoelastic case). Let u be the unique solution of (Q w ). (i) A unique p ∈ H 1 (ω) exists such that (u, p) satisfies (34), ( 35), (37)-( 39).

(ii) Moreover, if u ∈ H 1 (Ω), then (36) and (40) hold. In particular, (u, p) is the unique solution of (Q s ).

Proof. The result is stated using the same arguments that have been developped in the proof of Theorem 2.

In the next section, we provide some tools which enable the asymptotic equations of a viscoelastic flow in a thin domain to be solve. We present and analyse an algorithm and then, we focus on some applications which are related to lubrication theory. In particular, we illustrate boundary effects showing that the infinite journal bearing approximation, which is widely used in tribology, may lack relevance in viscoelastic regimes.

Numerical results and discusion.

As has been mentioned before, introducing viscoelastic effects leads to add a nonlinear term in the classical newtonian problem in pressure-velocity. This prevents us from following the classical method of obtaining only a problem in pressure both for the full continuous problem and for the numerical discretised one. So we have proposed a new method which will be presented in detail in section 3.3. This method is based on a two-step fixed point procedure. At the present time, we are not able to prove rigously the convergence of this method in the general case. However, we will give in section 3.2 a convergence result of each sub-step which can be considered as a new way of solving near-newtonian problem presented in section 3.1.

The method has been developed for a domain ω which is supposed to be rectangular with a size L × D.

3.1.

Numerical analysis for the Newtonian case. Let us recall the main equations of the Newtonian model:

(P)          - ∂ ∂z ∂ u ∂z + ∇ x p = 0, div x h 0 u(•, z) dz = 0.
In order to solve (P), a semi-discretized version of this problem, in the (x 1 , x 2 )-direction, is introduced. Thus, we use a centered structured grid based on a classical cell configuration (see Fig. 3). This particular case corresponds to an imposed flux on the left boundary x 1 = 0 and Dirichlet conditions for the pressure on the other boundaries. A similar discretisation may be adapted to the case of Dirichlet conditions for the pressure on the whole boundary. Let us denote by N = N x1 × N x2 the overall number of unknowns corresponding to this discretisation, by δ 1 (resp. δ 2 ) the step in the x 1 (resp.

x 2 ) direction, by h ij the value of h at a node (i, j). Furthermore, we denote

U (z) = (u ij (z)) i,j := (u(iδ 1 , jδ 2 , z)) i,j P = (p ij ) i,j := (p(iδ 1 , jδ 2 )) i,j
the semi-discretized horizontal velocity and discretized pressure.

Let A (resp. B) correspond to the x-discretisation of the operator ∇ (resp. div). We use the notation The problem (P) can be semi-discretised (i.e. discretised in the x-variable only) in the following way:

HU ij := hij 0 u ij (z) dz. u ij p ij v ij
(P ⋆ )      - ∂ ∂z ∂ U ∂z + A • P = 0, (42)
B • HU = 0. (43)
Concerning the boundary conditions, we impose that, for each node (i, j),

u ij ∈ H 1 (0, h ij ) with u ij (0) = s ij and u ij (h ij ) = 0,
the imposed velocity s being discetised by S = {s ij } ij . For the pressure, we impose a Dirichlet boundary condition which is written p ij = p ij for (i, j) at the boundary of the discrete domain, the imposed pressure is denoted by p ij . Notice that it is possible to solve (P ⋆ ) in a near analytic way by two integrations in the z-direction of the first equation in (P ⋆ ), taking into account the boundary condition on the velocity. We deduce

U = z (h -z) 2 A • P -S z -h h . (44)
Then putting the corresponding value of U as a function of the pressure P in the last equation of (P ⋆ ), we obtain the equation satisfied by P :

B • h 3 12 A • P = B • h 2 S . ( 45
)
This equation is the discretised finite difference formulation of the Reynolds equation whose solution P is unique and gives the knowledge of the velocity U by Equation (44).

As mentioned before, this last approach can not be generalised in the viscoelastic case. So we have proposed another algorithm which does not use the z-integration like the previous one. This algorithm is based on a fixed point formulation of the semidiscretised problem (P ⋆ ):

(P k )          - ∂ ∂z ∂ U k+1 (z) ∂z + A • P k = 0, (46)
P k+1 -P k + ρ B • HU k+1 = 0. ( 47
)
The stopping test of this process is based on the pressure error P k+1 -P k and on the velocity error U k+1 -U k . Note that the precision sought in pressure will induce a precision on the incompressibility condition via the parameter ρ. In fact, the algorithm is stopped as soon as |P k+1 -P k | is smaller than a prescribed value, noted r p (in the discrete ℓ 2 norm, for instance). If this condition is satisfied, it means in particular that the divergence term satisfies

max ij B • HU k+1 ij ≤ r p ρ ,
i.e. the free divergence equality is satisfied with an order of r p /ρ. For this reason r p /ρ will be called the "equilibrium parameter (for the free divergence condition)". In order to numerically attain the free divergence equality, we have to impose some value for r p satisfying r p ≪ ρ.

3.2.

Convergence of the method. We state the following theorem:

Theorem 5 (Convergence result). Assume that 0 < ρ < 3 h 3 L ∞ (ω) 1 δ 2 1 + 1 δ 2 2 .
Then for all k ∈ N, the problem (P k ) has a solution such that

(U k , P k ) ∈   ij H 1 (]0, h ij [)   × R Nx 1 ×Nx 2 .
Moreover a subsequence (still denoted {k}) exists such that, for all (i, j)

U k ⇀ U in ij H 1 (]0, h ij [), P k → P in R Nx 1 ×Nx 2 .
U and P being the solution of problem (P ⋆ ).

Proof. First of all, let us point out the fact that (P k ) is a linear problem. Thus, (U k , P k ) being given, (P k ) has a unique solution (U k+1 , P k+1 ). Using the linearity of problems (P k ) and (P ⋆ ), we prefer to work with the quantities U k = U k -U and P k = P k -P which satisfy problem (P k ) with homogeneous boundary Dirichlet conditions:

U k = 0, on Γ -∪ Γ + , (48)
P k = 0, on ∂ω. (49)
For the sake of simplicity in this proof, we write U k and P k instead of U k and P k (overscripts are dropped). At this point, we want to obtain estimates on the sequence (U k , P k ) and then prove that it converges to zero in appropriate spaces.

For each i, j, multiplying each component (46) ij by u k+1 ij , integrating over [0, h ij ] and then summing up for all i, j, we obtain, using an integration by parts,

∂ U k+1 ∂z , ∂ U k+1 ∂z ⋆ + A • P k , U k+1 ⋆ = 0, (50)
where (•,

•) ⋆ indicates the ⋆-scalar product (U , V ) ⋆ = i,j hij 0 u ij (z)v ij (z) dz
and, in the same way, (•, •) # is the scalar product defined by

(P, Q) # = i,j p ij q ij .
Now, | • | # being the associated norm, we deduce from Equation (47), after taking the #-scalar product by P k+1 , that

|P k+1 | 2 # -|P k | 2 # + |P k+1 -P k | 2 # + 2ρ B • HU k+1 , P k+1 # = 0. (51)
Moreover we obtain (first using a discrete integration by parts and then observing that P k+1 does not depend on z):

B • HU k+1 , P k+1 # = -A • P k+1 , HU k+1 # = -A • P k+1 , U k+1 ⋆ .
Using the previous equality, adding Equations ( 50) and (51) (with a multiplier 2ρ for (50)), we obtain the following estimate

|P k+1 | 2 # -|P k | 2 # + |P k+1 -P k | 2 # + 2ρ ∂ U k+1 ∂z 2 ⋆ = 2ρ A • (P k+1 -P k ), HU k+1 # . (52)
Now, we look for estimates for the right-hand side of this equality, denoted I 1 : As the operator A is bounded 1 by λ xy defined by

λ xy = 2 1 δ 2 1 + 1 δ 2 2 ,
we obtain that

I 1 ≤ 2ρ λ xy |P k+1 -P k | # HU k+1 # , 1
In fact, we obtain for instance

|A • P | 2 # = X ij " p i+1,j -p ij δ 1 « 2 + X ij " p i,j+1 -p ij δ 2 « 2 ≤ 2 " 1 δ 2 1 + 1 δ 2 2 « |P | 2 # .
and using the fact that for all (a, b) ∈ R 2 and α > 0, we obtain 2ab ≤ a 2 α + αb 2 , and we find that, for all α > 0

I 1 ≤ ρλ xy α |P k+1 -P k | 2 # + αρ HU k+1 2 # .
Moreover, for a regular function g : [0, h] → R such that g(h) = 0, we note that

g(z) = z h ∂ g ∂z (ξ) dξ so that, integrating over [0, h], h 0 g(z) dz = h 0 z h ∂ g ∂z (ξ) dξ dz = h 0 z ∂ g ∂z (z) dz,
and using the Cauchy-Schwarz inequality,

h 0 g(z) dz 2 ≤ h 0 z 2 dz h 0 ∂ g ∂z (z) 2 dz .
Thus we deduce that

HU k+1 2 # = ij hij 0 u k+1 ij (z) dz 2 ≤ C(h) ij hij 0 ∂ u k+1 ij ∂z (z) 2 dz,
where

C(h) = max i,j hij 0 z 2 dz = h 3 L ∞ (ω)

3

. That is,

HU k+1 2 # ≤ C(h) ∂ U k+1 ∂z 2 ⋆ .
We obtain

I 1 ≤ ρλ xy α |P k+1 -P k | 2 # + αρC(h) ∂ U k+1 ∂z 2 ⋆ . (53)
Putting Inequality (53) into (52), with an appropriate choice for α and ρ (to be detailed later), we can define two constants c 1 > 0 and c 2 > 0 such that:

|P k+1 | 2 # -|P k | 2 # + 1 - ρλ xy α c1>0 |P k+1 -P k | 2 # + ρ 2 -α C(h) c2>0 ∂ U k+1 ∂z 2 ⋆ ≤ 0. ( 54
)
The sign conditions on the two constants c 1 and c 2 are clearly satisfied if

0 < ρ < 2 C(h) λ xy := ρ crit. ,
and α being arbitrarily chosen in the set ]λ xy ρ, λ xy ρ crit. [. Notice that ρ crit. is a critical value of the parameter ρ allowing the above estimates. Summing the estimates (54) for k = 0 to k = K, we find bounds for U K and P K so that (U , P ) exists such that

U k ⇀ U in ij H 1 (]0, h ij [), P k → P in R Nx 1 ×Nx 2 .
Now, passing to the limit (K → +∞) in problem (P K ) enables us to deduce that (U , P ) is a solution of problem (P ⋆ ) (which has a unique solution, see Equation ( 45)). However, the boundary conditions, homogeneous (see Equations ( 48)-( 49)), imply that (U , P ) = (0, 0). This concludes the proof.

3.3. Numerical analysis for the viscoelastic case.

3.3.1.

Algorithm. The nonlinear problem (34)-( 40) is solved using a fixed-point method at different levels of the resolution. Let us define a continuous fixed-point procedure. The idea of the general algorithm relies on the possibility of reaching a solution to the nonlinear problem (34)-(40) as the limit (n → +∞) of the following problem:

(P n )          - ∂ ∂z f (u n ) ∂ u n+1 ∂z + ∇ x p n+1 = 0, div x h 0 u n+1 (•, z) dz = 0. with f (u n ) = (1 -r) + r 1 + De 2 (1 -a 2 ) ∂ u n ∂z 2
In order to solve (P n ), the same semi-discretisation in the (x 1 , x 2 )-direction as in the newtonian case, is used. Now, we present the algorithm which solves the semi-discrete version of (P n ). The way of computing U n+1 and P n+1 is provided by the algorithm presented in the Newtonian case :

                       Input: U n,0 = U n , P n,0 = P n , Loops on k: (P k n )        - ∂ ∂z f (U n ) ∂ U n,k+1 ∂z + A • P n,k = G, P n,k+1 -P n,k + ρ B • HU n,k+1 = 0, Output: U n+1 = U n,∞ , P n+1 = P n,∞ .
The algorithm is stopped as soon as P n+1 -P n is smaller than a prescribed value (in the discrete ℓ 2 norm, for instance).

3.3.2.

Remarks on the method. The algorithm that we propose views the viscoelastic problem as a sequence of Newtonian-type problems. Formally, the numerical solution attained is a fixed-point solution of the semi-discretized version of (P n ).

Following the same idea as in the Newtonian case, the theoretical study establishes the boundedness of the sequence, provided some constraints (which do not depend on k and n) are respected. More precisely we can notice that, since the function f satisfies f ≥ (1-r), then we obtain estimates which do not depend on n and k under the condition

0 < ρ < 3 (1 -r) h 3 L ∞ (ω) 1 δ 2 1 + 1 δ 2 2 .
This condition is more restrictive than the preceeding one but sufficient for all the n-steps.

Unfortunately it is not so obvious that the sequence of solutions (U n , P n ) converges to a fixed-point solution of the semi-discretised version of problem P n , because of the nonlinearity which leads to a lack of compactness. However, in practical situations, we observe the following phenomena: (i) Under the constraint r < 8/9, we observe that the algorithm converges to a numerical viscoelastic solution under the above condition. (ii) Under the constraint r > 8/9, it is observed that the sequence of the numerical Newtonian-type solutions does not converge to a viscoelastic one, which may be related to the non-uniqueness result (see [START_REF] Boyer | Numerical study of viscoelastic mixtures through a Cahn-Hilliard flow model[END_REF] for similar observations in a Stokes-Oldroyd flow).

3.4. Numerical results. In this subsection, we propose three series of numerical tests:

• Test 1: we study the influence of the numerical parameters on the solution. In particular, the control of ρ with respect to the stopping error may have some influence on the numerical solution. However, at least in the Newtonian case, we illustrate the behaviour of the solution with respect to ρ and show that it converges to the solution of the Reynolds newtonian equation (obtained by solving Equation ( 45)) as ρ tends to 0. This solution will be denoted "Reynolds" in the following.

• Test 2: we study the influence of the Deborah number.

• Test 3: we show that three-dimensional effects may occur. In particular, the approximation of the "journal bearing of infinite width", which is valid (and widely used) in the Newtonian case, cannot be considered due to viscoelastic effects. For this, the following data has been used (see Table 1).

3.4.1. Test 1: influence of the numerical parameters. In this setting, we study the purely Newtonian case, which allows us to compare our numerical pressure to the theoretical one: the solution of the classical Reynolds equation. In particular, we focus on the role of the equilibrium parameter r p /ρ (corresponding to the error on the free divergence condition). Since our goal is to obtain simultaneously the convergence of the pressure Test 1

Test 2 Test 3

Domain ω and the equilibrium of the free divergence condition, we first impose an artificial time step ρ = 10 -3 , which ensures the convergence of the method. Then, we choose different values of r p /ρ in order to observe its numerical influence over the corresponding solution.

[0, 1] × [0, 5] [0, 1] × [0, 5] [0, 1] × [0, 5] Gap h(x) (2x 1 -1) 2 + 0.5 1 -0.3x 1 + 0.5x 2
In particular, it is sufficient to compare our numerical solution (for different values of r p /ρ) to the Reynolds one. Thus we consider numerical data given in Table 1. Let us state the values for the boundary conditions; at the (left) boundary x 1 = 0, the normalized flux is given by q 0 = 0.3s x h |x1=0 while, at other boundaries, the pressure is p = 0. Now the influence of the ratio r p /ρ is illustrated in Fig. 4; the left-hand side figure is the Reynolds pressure distribution, in the full domain. The right-hand side figure enables in the x 1 direction (at a fixed x 2 , namely x 2 = x 0 2 = 2.5) the solutions to be observed corresponding to different values of r p /ρ. It can be observed that the numerical pressure tends to the Reynolds one as the value of r p /ρ decreases. At r p /ρ = 10 -4 , numerical and Reynolds solutions even coincide. To be mentioned is the peculiar slope near the inlet to the specific data for x 1 = 0 which is related to the input flux and not to the pressure as on the other part of the boundary. In this subsection, we compare our model to the ones developed by F.T. Akyildiz and H. Bellout [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF] and J.A. Tichy [START_REF] Tichy | Non-Newtonian lubrication with the convective Maxwell model[END_REF]. Notice that, unlike our model, these previous works only deal with two-dimensional flows, corresponding, for example, to journal bearings with an infinite width (i.e. devices whose size satisfy D/L > 4). This assumption enables that, up to boundary effects localized at x 2 = 0 and x 2 = D, to consider that the flow is mainly described by its behaviour at a cross section (x 2 = D/2 for instance) and that it remains the same at another cross sections (as long as it is far from the boundaries). Following the work of F.T. Akyildiz and H. Bellout [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF], we choose the physical data given at Table 1. To complete the scope of the boundary conditions, let us mention that p = 0 is imposed on the whole boundary ∂ω.

More precisely, in order to observe the effects of the Deborah number over the pressure distribution, we have used the same values as in the paper of F.T. Akyildiz and H. Bellout [START_REF] Talay Akyildiz | Viscoelastic lubrification with Phan-Thein-Tanner fluid (PTT)[END_REF]: De = 0.1, De = 0.2,..., De = 3. We can observe the behaviour of the solution, as De increases, in Fig. 5, corresponding to the pressure profiles at a fixed x 2 = 2.5.

For small values of De (0, 0.1, 0.2, 0.3), the results are similar to the ones of F.T. Akyildiz and H. Bellout, except that they have been generalized to a three-dimensional flow: viscoelastic effects tend to dampen the peak pressure.

For large values of De (1, 2,...), the results differ from the ones of F.T. Akyildiz and H. Bellout. This is due to the fact that our initial models are different (and so are the corresponding asymptotic analyses). Here, we may observe that the viscoelastic solution, as De → +∞, converges to the solution of the purely viscous solution with an relative viscosity parameter (1r) (instead of 1). The viscoelastic nonlinear contribution formally tends to vanish for large values of De. However our model is not relevant for large values of De, since De is assumed to be of the order of ε. As was pointed out in the previous subsection, when the length and width of a device satisfies D/L > 4, a classical approximation is used in lubrication theory. The one-dimensional Reynolds equation is used to describe the behaviour of the flow at any cross-section which is not located at boundaries x 2 = 0 or x 2 = D. This assumption enables the space dimension in the analysis of such phenomena to be reduced. This is well understood in the Newtonian case but numerical tests illustrate that such an assumption is not necessarily relevant when viscoelastic effects occur. In fact, three dimensional boundary layers are induced by viscoelastic effects.

We have used the physical and numerical data given in Table 1. Let us give the values for the boundary conditions; at the (left) boundary x 1 = 0, the normalized flux is given by q 0 while, at other boundaries, the pressure is p = 0.

In Fig. 6, (from left to right, top to bottom), we have the pressure profiles corresponding to r = 0 (Newtonian case), r = 0.2, r = 0.5 and r = 0.8. Obviously the one-dimensional flow assumption in the Newtonian case is valid as long as the crosssection is not located at the boundaries but we can see that this assumption does not hold anymore for increasing values of r.

Fig. 2 .

 2 Fig. 2. Mixed boundary conditions

  and we verify that ϕ ∈ K(0, 0) ∩ D(Ω) and Φ(ϕ) = f . ⊲ Step 2-2. (Constant pressure on each connected component of ∂ω p ) -Let us defineZ = {ϕ ∈ D(∂ω) such that ∂ϕ ∂τ = 0 on ∂ω q }(τ being the tangent vector to the boundary ∂ω). For all ϕ ∈ Z, we extend ϕ on ω and define f = rot ϕ. Since div(rot) = 0 and n• rot = ∂/∂τ , we deduce that f ∈ X. So ψ ∈ K(0, 0) ∩ D(Ω) exists such that rot ϕ = h 0 ψ(•, z) dz.From (31), we deduceω ∇ x p • rot ϕ = 0.After integrating by parts, we obtain ∀ϕ ∈ Z,

Fig. 3 .

 3 Fig. 3. Spatial discretisation and position of the unknowns

Fig. 4 .

 4 Fig. 4. Influence of the equilibrium parameter

Fig. 5 .

 5 Fig. 5. Influence of the Deborah number

Fig. 6 .

 6 Fig. 6. Pressure profiles for r = 0.0, r = 0.2, r = 0.5 and r = 0.8

Table 1 .

 1 Numerical data

	1 (2x 1 -1) 2 + 0.5