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

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

Carlos Vázquez
email: carlosv@udc.es

Micro-Roughness Effects in (Elasto)Hydrodynamic Lubrication Including a Mass-Flow Preserving Cavitation Model

Keywords: Elastohydrodynamic, Lubrication, Roughness

An average Reynolds equation is proposed for predicting the effects of deterministic periodic roughness, taking JFO mass flow preserving cavitation model and elastohydrodynamic effects into account. For this, the asymptotic model is based upon a double scale analysis approach. The average Reynolds equation can be used both for the description of cavitation on a macroscopic scale in widening gap regions as well as for cavitation on the microscopic interasperity scale. Results of numerical simulations based on the model are presented for the case of a hydrodynamically lubricated journal bearings and an elasto-hydrodynamically lubricated point contact.

In this paper, it is explained how the double scale procedure, already used to obtain average equations with periodic roughness in the case of rigid bearings [START_REF] Bayada | A double-scale analysis approach of the Reynolds roughness. Comments and application to the journal bearing[END_REF][START_REF] Bayada | Effets d'anisotropie par homogénéisation dans un problème à frontière libre[END_REF][START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF], can be extended to EHD problems including cavitation and starvation.

The JFO mass flow preserving model is used, including pressure and saturation as unknown functions. This model takes into account both microcavitation (due to the microroughness) and macrocavitation (due to the diverging part of the gap). The average equation can be easily solved for some specific roughness patterns (transverse, longitudinal) exactly in the same way as the initial EHD problem with cavitation. Numerical results are given for both purely hydrodynamic and EHD point-contact problems, for a two dimensional device.

Basic equations

Our studied cavitation model, like the Elrod algorithm and its variants [START_REF] Kumar | A finite element cavitation algorithm[END_REF][START_REF] Shi | A mixed soft elastohydrodynamic lubrication model with interasperity cavitation and surface shear deformation[END_REF], views the film as a mixture. It does not, however, make the assumption of liquid compressibility in the full film area as in [START_REF] Vijayaraghavan | An efficient, robust, and time accurate numerical scheme applied to a cavitation algorithm[END_REF] and some other papers [START_REF] Elrod | A cavitation algorithm[END_REF].

In this general framework, the effective gap between the two close surfaces, denoted h [p] , contains a given rigid contribution h r and an elastic one, which strongly depends on the main unknown p (lubricant pressure) in the following nonlocal form:

h [p] (x) = h r (x) + Ω k(x, z)p(z) dz,
the kernel k depending on the kind of contact. The classical approximation of the rigid gap [START_REF] Dowson | Elastohydrodynamic lubrication[END_REF] is given by the expression

h r (x) =        h 0 + x 2 1 + x 2 2

2R

, for point contacts,

h 0 + x 2 1 R , for line contacts, (1)
that represents a parabolic approximation for a given sphere-plane (point contact) or cylinder-plane (line contact) gap, R being the sphere or cylinder section radius. In the point contact case, the positive constant h 0 corresponds to the gap at the point nearest to contact while, in the line contact case, it includes a constant related to the logarithmic kernel (see [START_REF] Johnson | Contact mechanics[END_REF] for instance). Now, let us introduce the general property of the kernel k:

k(x, z) =          c 0 log c 1 -z 1 x 1 -z 1 , for line contacts, c 0 (x 1 -z 1 ) 2 + (x 2 -z 2 ) 2
, for point contacts,

where c 0 > 0 and c 1 ≥ max{ x 1 , x ∈ Ω}. The flow obeys the following "universal" Reynolds equation (here written in a dimensionless form):

2 i=1 ∂ ∂x i h 3 [p] e -αp ∂p ∂x i = ∂ θh [p] ∂x 1 , (3)
p ≥ 0, 0 ≤ θ ≤ 1, p (1 -θ) = 0, (4)
p is the pressure (assumed to be a positive function), θ is the relative mixture density, h [p] the real film thickness, x 1 is the direction of the effective relative shear velocity of the device, while x 2 is the transverse direction. Here, the pressure cavitation is assumed to be zero. Moreover, the lubricant is piezoviscous so that the viscosity obeys the Barus law (notice that other laws may be taken into account), α being the piezoviscosity parameter.

3 Asymptotic expansion

General case in the elastohydrodynamic configuration

Let us suppose that the roughness is periodically reproduced in the two x 1 and x 2 directions from an elementary cell Y (or "miniature bearing" in Tonder's terminology). Relatively to this elementary cell Y , we define the microscopic

εY Y y 1 y 2 x 1 x 2 y = x ε
Fig. 1. Macroscopic domain and elementary cells variables y = (y 1 , y 2 ) such that Y = {0 < y 1 < 1, 0 < y 2 < y 2 } which is used to describe the oscillatory effects (see Fig. 1).

Let us now consider gaps that can be written as h r (x, y). In this section, we use the two-scale analysis which has been widely developped for the hydrodynamic case in [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF]. The equations should be read as

2 j=1 ∂ ∂x j h [p] (x, y) 3 e -αp(x) ∂p ∂x j (x) = ∂ ∂x 1 θ(x) h [p] (x, y) , p(x) ≥ 0, 0 ≤ θ(x) ≤ 1, p(x) (1 -θ(x)) = 0,
with the additional notation:

h [p] (x, y) = h r (x, y) + Ω k(x, z) p(z) dz.
Let us introduce the parameter ε which denotes the period of the roughness patterns: the fast variables (y 1 , y 2 ) describe oscillatory effects induced by the roughness patterns, namely

(y 1 , y 2 ) := x 1 ε , x 2 ε
Thus, similarly to the hydrodynamic case (studied in [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF]), we shall look for an asymptotic expansion of the solutions:

p(x) = p 0 (x) + εp 1 (x, y) + ε 2 p 2 (x, y) + ..., (5) θ
(x) = θ 0 (x, y) , (6)
p 0 being a positive function, and each unknown p i (i ≥ 1) and θ 0 being function of (x, y). The problem of the boundary conditions to be satisfied by the p i is somewhat difficult but may be summarized so:

(i) The natural boundary conditions on (p, θ) are assigned to p 0 and an equivalent saturation linked to θ 0 . (ii) The functions p i , i ≥ 1, are Y periodic, i.e. periodic in the two variables y 1 , y 2 , for each value of (x 1 , x 2 ).

To be noticed that, unlike of p, we do not introduce an asymptotic expansion for θ. This can be explained by observing the evolution of p and θ as ε tends to 0. Clearly, the oscillations of the pressure are decreasing and p tends to a smooth function (namely p 0 which, actually, does not depend on the fast variable). This is not the case for θ and an asymptotic smooth limit cannot be considered. Moreover, from the initial properties of the pressure-saturation, the following ones hold:

p 0 (x) ≥ 0, 0 ≤ θ 0 (x, y) ≤ 1, p 0 (x, y) (1 -θ 0 (x, y)) = 0. ( 7
)
Using the two-scale analysis [START_REF] Bayada | A double-scale analysis approach of the Reynolds roughness. Comments and application to the journal bearing[END_REF][START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF], we obtain the micro-macroscopic link between the unknowns p 0 , p 1 and θ 0 :

p 1 (x, y) = e αp 0 (x) χ (1,0) [p 0 ] (x, y) - 2 i=1 ∂p 0 ∂x i (x) w (i) [p 0 ] (x, y), (8)
in which w

(i) [p 0 ] , χ (i),0 [p 0 ]
(and, additionally, χ

(i,⋆) [p 0 ] ) (i = 1, 2
) are the Y periodic solutions (up to an additive constant) of the following local problems (see [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF] for the rigorous introduction of these local problems):

2 j=1 ∂ ∂y j   h 3 [p 0 ] e -αp 0 ∂w (i) [p 0 ] ∂y j   = ∂h 3 [p 0 ] ∂y i , (9)
2 j=1 ∂ ∂y j   h 3 [p 0 ] e -αp 0 ∂χ (i,0) [p 0 ] ∂y j   = ∂θ 0 h [p 0 ] ∂y i , ( 10
)
2 j=1 ∂ ∂y j   h 3 [p 0 ] e -αp 0 ∂χ (i,⋆) [p 0 ] ∂y j   = ∂h [p 0 ] ∂y i . ( 11
)
Notice that, contrary to the local problems in the purely hydrodynamic case which have been studied in [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF], the local problems in the EHL problem highly depend on the macroscopic pressure p 0 , which is due to the piezoviscosity and elastic deformation. These local problems describe the link between the unknowns at the microscopic scale, the macroscopic one playing the role of a parameter. To complete the results, we obtain the following average Reynolds equation, taking into account the roughness effects. There are now three leading unknowns (p 0 , Θ 1 , Θ 2 ) satisfying the following equations:

2 i,j=1 ∂ ∂x i A ⋆ [p 0 ] ij e -αp 0 ∂p 0 ∂x j = 2 i=1 ∂ ∂x i Θ i B (i,⋆) [p 0 ] , (12)
p 0 ≥ 0, p 0 (1 -Θ i ) = 0, (i = 1, 2), (13)
where the average coefficients are given by:

A ⋆ [p 0 ] ij = h 3 [p 0 ] δ ij -h 3 [p 0 ] ∂w (i) [p 0 ] ∂y j Y , B (i,⋆) [p 0 ] = h [p 0 ] δ 1i -h 3 [p 0 ] ∂χ (1,⋆) [p 0 ] ∂y i Y ,
δ ij denoting the Kronecker symbol (δ ij has value 1 for i = j, 0 for i = j). Equations ( 12) and ( 13) deal with any periodic roughness pattern. To be noticed is the fact that the differential operator is no more of the Reynolds type since extra terms ∂ 2 p 0 /∂x i ∂x j appear. The right-hand side also contains an additive term in the x 2 direction. Moreover, it is possible to prove in a rigorous way [START_REF] Bayada | Homogenization of a nonlocal elastohydrodynamic lubrication problem: a new free boundary model[END_REF] that there exists a solution with an isotropic saturation Θ = Θ 1 = Θ 2 and 0 ≤ Θ ≤ 1. Under some additional assumptions (see [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF]), we can even prove some complementary results, which are presented in the next subsections.

Particular cases in the elastohydrodynamic configuration

When dealing with transverse roughness and longitudinal roughness, the average coefficients can be easily simplified and given in an explicit form. In both cases, the asymptotic system has the same structure than the initial (smooth) one, i.e.

2 i=1 ∂ ∂x i A (i,⋆) [p 0 ] e -αp 0 ∂p 0 ∂x i = ∂ΘB (1,⋆) [p 0 ] ∂x 1 , (14)
p 0 ≥ 0, 0 ≤ Θ ≤ 1, p 0 (1 -Θ) = 0. (15)
⊲ Transverse roughness: the average coefficients are

A (1,⋆) [p 0 ] = 1 h -3 [p 0 ] Y , A (2,⋆) [p 0 ] = h 3 [p 0 ] Y , B (1,⋆) [p 0 ] = h -2 [p 0 ] Y h -3 [p 0 ] Y ,
the link between the macroscopic saturation Θ and the microscopic one θ 0 being given by:

Θ = 1 h -2 [p 0 ] Y θ 0 h 2 [p 0 ] Y .
⊲ Longitudinal roughness: the average coefficients are

A (1,⋆) [p 0 ] = h 3 [p 0 ] Y , A (2,⋆) [p 0 ] = 1 h -3 [p 0 ] Y , B (1,⋆) [p 0 ] = h [p 0 ] Y ,
the link between the macroscopic saturation Θ and the microscopic one θ 0 being given by:

Θ = θ 0 h [p 0 ] Y h [p 0 ] Y .
All the earlier results are valid for both elastohydrodynamic and hydrodynamic cases and, thus, generalize the ones that have been stated in [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF]. As an important feature, Θ is not the average of the microscopic saturation θ 0 but contains some anisotropic effects due to the roughness direction.

Particular cases in the hydrodynamic configuration

In the purely hydrodynamic case, one can prove some additional results, corresponding to a wide class of two dimensional roughness patterns. Indeed, suppose that h r can be written under the form

h r (x, y) = h 1 (x, y 1 ) h 2 (x, y 2 ) ,
then we get the following (hydrodynamic) average equation

2 i=1 ∂ ∂x i A (⋆) ii ∂p 0 ∂x i = ∂ ∂x 1 ΘB (1,⋆) , p 0 ≥ 0, 0 ≤ Θ ≤ 1, p 0 (1 -Θ) = 0, with A (⋆) ii = h 3 j Y h -3 i Y , for i, j = 1, 2 and i = j, B (1,⋆) = h -2 1 Y h -3 1 Y h 2 Y ,
the link between the micro-macroscopic saturations being given by

Θ = 1 h 2 Y h -2 1 Y θ 0 h 2 h 2 1 Y . 8
In this section, the numerical simulation of a micro(elasto)hydrodynamic contact is performed to illustrate the theoretical results of convergence stated in the previous sections. For the numerical solution of the ε-dependent problems and their corresponding average one, the algorithm is based on the Bermudez-Moreno method : a characteristics method coupled to a duality method, adapted to steady state problems dealing with a finite element spatial discretization, is used in order to treat the nonlinearity related to cavitation phenomena. Piezoviscous effects and elastic deformation of the surfaces are taken into account by a standard fixed-point method. These numerical techniques have been already successfully used in previous papers dealing with hydrodynamic aspects (see [START_REF] Bayada | Characteristics method for the formulation and computation of a free boundary cavitation problem[END_REF][START_REF] Bermúdez | La méthode des caractéristiques pour les problèmes de convection-diffusion stationnaires[END_REF]), and elastohydrodynamic aspects (see, for instance, [START_REF] Arregui | Mathematical analysis and numerical simulation of a Reynolds-Koiter model for the elastohydrodynamic journalbearing device[END_REF][START_REF] Durany | Numerical simulation of a lubricated Hertzian contact problem under imposed load[END_REF]). It has been proved that, for small loads, it is a rigorous method whose numerical solution converges to the continuous one, as the mesh size tends to 0. In Table 1, we present the functional coefficients

A (1,⋆) [p 0 ] , A (2,⋆) [p 0 ] , B (1,⋆) [p 0 ]
that appear in the average problem for purely transverse and purely longitudinal roughness cases which have been partially computed with MAPLE. Suppose that the domain to be considered is a rectangular one (a, b) × (c, d), and that the rough (elasto)hydrodynamic gap can be written under the form

h [p] x, x ε = h r (x)+h d [p](x)+        α 1 sin 2π ε x 1 -a b -a , (transverse roughness) α 2 sin 2π ε x 2 -c d -c , (longitudinal roughness)
where h r (resp. h d [p]) denotes the smooth rigid (resp. elastic) contribution to the gap and the remaining term describes the roughness patterns in the purely transverse or longitudinal case.

Transverse roughness

Longitudinal roughness

h [p] (x, y) h r (x) + h d [p](x) + α 1 sin (2πy 1 ) h r (x) + h d [p](x) + α 2 sin (2πy 2 ) A (1,⋆) [p 0 ] 2 (h r + h d [p]) 2 -α 2 1 5/2 2(h r + h d [p]) 2 + h 2 r (h r + h d [p]) 3 + 3 2 (h r + h d [p]) α 2 2 A (2,⋆) [p 0 ] (h r + h d [p]) 3 + 3 2 (h r + h d [p]) α 2 1 2 (h r + h d [p]) 2 -α 2 2 5/2 2(h r + h d [p]) 2 + α 2 2 B (1,⋆) [p 0 ] 2(h r + h d [p]) (h r + h d [p]) 2 -α 2 1 2(h r + h d [p]) 2 + α 2 1 h r + h d [p]
Table 1 Elastohydrodynamic average coefficients

Hydrodynamic case

We adress the numerical simulation of journal bearing devices with circumferential supply of lubricant. The mechanical characteristics of the device are given by:

• length: L = 0.019 m,

• bearing radius: R b = 0.0164975 m,

• journal radius: R j = 0.01647 m,

• mean radius: R m = 0.5 (R b + R j ),

• clearance : c = R b -R j ,
• eccentricity:

ρ = 0.2.
The physical characteristics of the regime are the following ones:

• supply pressure: p a = 283000 P a,

• lubricant viscosity: µ = 0.02 P a.s,

• shear velocity: v 0 = 17.247 m/s.

The earlier problem leads to the following set of equations

2 i=1 ∂ ∂x i h 3 ∂p ∂x i = Λ ∂ θh ∂x 1 , p ≥ 0, 0 ≤ θ ≤ 1, p (1 -θ) = 0,
with Λ = 6µv 0 , and the real smooth gap shoud be read as

h(x) = c (1 + ρ cos (x 1 /R m )) .
The equations have to be solved on the domain (0, 2πR m ) × (0, L/2) (see the developped configuration on Fig. 2) with the following boundary conditions:

• p = p a on Γ a , p = 0 on Γ 0 ,

• conditions on Γ ♯ : p and Λ θh -h 3 ∂p ∂x 1 are 2πR m x 1 periodic.

To be noticed is the fact that the corresponding average boundary conditions can be written as follows (using the terminology of Table 1):

• p 0 = p a on Γ a , p 0 = 0 on Γ 0 ,

• conditions on Γ ♯ : p 0 and Λ ΘB (1,⋆) -A (1,⋆) ∂p 0 ∂x 1 are 2πR m x 1 periodic. Actually, we consider transverse roughness patterns and the gap should be read as:

x 1 x 2 Γ ♯ Γ ♯ Γ a Γ 0 2πR m L 2
h x, x ε = c 1 + ρ cos (x 1 /R m ) + a r sin 2π ε x 1 R m
with a r /ρ = 0.2, a r denoting the amplitude parameter of the roughness. Numerical tests have been made with a mesh grid containing 201*51 elements, which is proved to be sufficient in terms of numerical convergence. As there is at most 30 roughness patterns in the x 1 direction, this has lead to a mesh discretization of 7 points at least for each roughness pattern. Moreover, numerical results remain the same when using meshes with size 401*51 or 201*101.

Fig. 3 and 4 (resp. 5 and 6) show the pressure and saturation profiles for ε = 1/15 (resp. ε = 1/30) compared to the average solution, at a fixed x 0 2 = L/4. Thus, it allows to observe the roughness effects in the x 1 direction. The amplitude of the pressure oscillations tend to be damped, although the amplitude of the saturation oscillations stay the same in cavitated areas. As it was noticed in [START_REF] Bayada | An average flow model of the Reynolds roughness including a mass-flow preserving cavitation model[END_REF], it points out the fact that when the number of roughness patterns increases, the pressure behaves as a smooth function, namely p 0 (x), while the saturation behaves as a highly oscillating function, namely θ 0 (x, x/ε). Thus, the pressure tends to a smooth one as ε tends to 0, while the saturation is always oscillating. To be noticed on Fig. 4 and6 is the fact that the cavitation area is made of two macrocavitation zones (for ε = 1/15, x 1 > 0.06 and x 1 < 0.01) and a lot of microcavitation zones. Fig. 7 and8 represent the average pressure and saturation in the real domain.

Let us notice that the average process, that we have presented here, only allows a macroscopic description of the microscopic roughness effects as the number of patterns is large: in particular, interasperity cavitation cannot be predicted in the average model, which only describes the (macro-) cavitation regions for a large number of roughnesses. Additional information on interasperity cavitation is contained in the supplementary terms of the asymptotic expansion of the pressure, namely p 1 (x, y), for which little information is known.

Elastohydrodyamic case

The numerical tests deal with a dimensionless problem, as described in Section 2: the domain is (-4, 2) × (2, 2). The considered rigid contribution to the gap is a normalized one:

h 0 + x 2 1 + x 2 2
2 where h 0 denotes the minimum thickness. Since a point contact has been considered, we choose the following Hertz model:

k(x, z) = 2 π 2 1 (x 1 -z 1 ) 2 + (x 2 -z 2 ) 2
The equations are:

2 i=1 ∂ ∂x i h[p] 3 e -αp ∂p ∂x i = ∂ θh[p] ∂x 1 , p ≥ 0, p (1 -θ) = 0 0 ≤ θ ≤ 1,
The boundary conditions are the following ones:

x 1 x 2 Γ ⋆ 6 4
Fig. 9. Normalized EHL domain

• flow condition on Γ ⋆ : θh[p] -h[p] 3 ∂p ∂x 1 = θ in h[p], with θ in = 0.3,
• p = 0 elsewhere.

The corresponding average boundary conditions can be written as:

• flow condition on Γ ⋆ : ΘB

(1,⋆) [p 0 ] -A (1,⋆) [p 0 ] ∂p 0 ∂x 1 = θ in C (⋆) [p 0 ] , with C (⋆) [p 0 ] = (h r + h d [p]) |Γ⋆ . • p 0 = 0 elsewhere.
The chosen values of h 0 and α will be discussed further.

Transverse roughness

Numerical tests deal with the following rigid contribution to the gap:

h 0 + x 2 1 + x 2 2 2 + h 0 sin 2π x 1 + 4 6 ε
with h 0 = 0.5 and different values of ε. Moreover, the piezoviscosity has been taken to α = 1. This numerical configuration corresponds to a very small load contact. Computations have been made for different values of ε (with at most 1/ε = 30 roughness patterns), and the mesh grid which has been used contains 201*51 elements, as in the hydrodynamic case. The average coefficients in the transverse roughness cases are deduced from Table 1 (see page 9).

Fig. 10-11 (resp. 12-13) provide the pressure-saturation profiles at x 0 2 = 0 with ε = 1/20 (resp. ε = 1/30) compared to the average solution. To be observed is the fact that the average profiles give a satisfying approach of the roughness effects (ε = 1/30): indeed, the pressure profiles given in Fig. 10 and 12 evidence the fact that the average pressure is a smooth version of the rough pressure. Similarly, by Fig. 11 and 13, the average saturation can be seen as an average version of the rough saturation, up to anisotropic effects. On Fig. 16, we observe that the deformation corresponding to ε = 1/30 nearly coincides with the average one: this is due to the regularizing effects of the Hertz kernel. In fact, the deformation profile has a rate of convergence which is much greater than the pressure profile. Fig. 14, 15 and 17 represent the average pressure, saturation and deformation in the domain.

Again, similarly to the hydrodynamic computations, the average model does not provide any information of cavitation occurring at the microscopic level: it only means that the effective pressure is closer and closer to the average (smooth) pressure as the number of roughness patterns increases (i.e. ε tends to 0) but it does not provide information of micro-cavitation for a given value of ε.

Longitudinal roughness

Numerical tests have been made for the following rigid contribution to the gap:

h 0 + x 2 1 + x 2 2 2
+ h 0 sin 2π x 2 + 2 4 ε with h 0 = 0.5 and different values of ε. Moreover, the piezoviscosity has been taken to α = 1. Let us notice that, again, average coefficients in the transverse roughness cases are deduced from Table 1 (see page 9).

Computations have been made for different values of ε (with at most 1/ε = 30 roughness patterns), with a mesh grid which containing 51*201 elements: the discretization has been refined in the x 2 direction in order to take into account the roughness effects in this direction.

Fig. 18 and 19 represent the pressure and deformation profiles at x 0 1 = -0.4, in the x 2 direction (in order to observe the roughness effects). This choice corresponds to the maximum pressure in the average case, which is attained at (x 0 1 , x 0 2 ) = (-0.4, 0). Of course, the saturation profile is omitted, for all corresponding saturation functions would be identically equal to 1 (no cavitation in this part of the domain). Significantly, the size of the oscillations for the pressure are damped easily, and convergence of the rough solution to the average one is illustrated on both figures. Similarly, pressure / saturation / deformation curves are omitted, for they are similar to the ones observed in the transverse roughness case.

Influence of the roughness effects in EHL and hydrodynamic cases

Fig. 20 and 21 show the difference of roughness effects between a purely hydrodynamic (isoviscous) configuration and an elastohydrodynamic (piezoviscous) configuration. The numerical and physical data are the same as in the transverse roughness case (see Subsection 4.2.1), except for the rough gap whose amplitude of roughness patterns is modified in order to prevent contact between the surfaces in the hydrodynamic case: thus, the gap is

h 0 + x 2 1 + x 2 2 2 + 0.7 h 0 sin 2π x 1 + 4 6 ε ,
all other numerical data being the same as before (in particular for the value of h 0 and the piezoviscosity parameter α). It can be noticed that the elevation of the pressure due to the roughness patterns is less important in the EHL case than in the purely hydrodynamic case. This is due to the fact that the elastic deformation tends to damp the additional load corresponding to the roughness. It has little influence over the saturation distribution, although the homogenization process does not allow to get microcavitation effects which do exist when a deterministic rough pattern is considered. Though, this analysis also states that microcavitation effects tend to vanish as ε tends to 0.

Influence of the roughness over the load

Numerical tests have been made for the following rigid contribution to the gap:

h 0 + x 2 1 + x 2 2 2
+ a r sin 2π

x 1 + 4 6 ε with h 0 ∈ {0.5, 1, 1.5, 2} and a r /h 0 ∈ {0.2, 0.4, 0.6, 0.8, 1}. Moreover, the elastic contribution to the gap is the one given at the beginning of Subsection 4.2, and piezoviscosity has been taken to α = 0 (isoviscous case). Results have been obtained with a mesh grid 101*101. They are given on Fig. 22, showing the influence of the minimum thickness h r -a r over the load W for different values of h 0 . Results are taken from the analysis of the corresponding average solution.

Influence of the piezoviscosity

We focus on the behaviour of the solution with respect to the piezoviscosity parameter α. Numerical and physical data are the same than in the elastohydrodynamic case with transverse roughness (Subsection 4.2.1), except that we take into account piezoviscous properties of the lubricant: α = 0, 1, 2 or 3.

Fig. 23 and 24 represent the pressure and deformation profiles at x 0 2 = 0 in the average case. They illustrate the trend induced by the piezoviscosity parameter: the peak pressure and the peak deformation increase with α. Only a few variations affect the saturation distribution (for this reason, the corresponding curves are omitted).

Table 2 (resp. Table 3) gives the variation of the peak pressure (resp. deformation) with respect to the isoviscous case (α = 0) in different rough cases (including smooth and average ones). The relative variation of the peak pressure is denoted ∆p/p = max(p α ) -max(p 0 ) max(p 0 ) where p α (resp. p 0 ) denotes the pressure distribution corresponding to the piezoviscous regime α = 0 (resp. isoviscous regime α = 0). Similarly, the 3 Maximum deformation elevation due to piezoviscosity

Conclusion

Numerical simulations of EHL problems are important for the design of both components and lubricants. On the one hand, th key numerical methods used in EHL solvers are the multigrid method [START_REF] Venner | Multi-level methods in lubrication[END_REF] and multilevel multi-integration [START_REF] Brandt | Multilevel matrix multiplication and fast solution of integral equations[END_REF], which have proved robust and highly desirable in terms of reducing the time spent in solving EHL problems, at least for a moderated number of roughness patterns. On the other hand, accurate solutions of micro-EHL problems with realistic surface roughness requires such a high spatial resolution that even a single solution may be a large computational task [START_REF] Lubrecht | Elastohydrodynamic lubrication of rough surfaces[END_REF]. Thus, as ε tends to 0, it seems more relevant to use the solution procedure that we have proposed because it deals with smooth coefficients. However, this average method may be coupled to multigrid methods (especially for the computation of the coefficients at the lowest grid resolution). The average method is valid for transverse or longitudinal roughness patterns. Further investigation has to be made in order to take into account anisotropic two dimensional effects.

  with respect to y L = length of the journal bearing R b , R j , R m = bearing radius, journal radius, mean radius Γ a , Γ 0 , Γ ♯ , Γ ⋆ = boundaries of the device c, ρ = clearance, eccentricity of the bearing

Fig. 2 .

 2 Fig. 2. Journal bearing domain

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .

 345678 Fig. 3. Hydrodynamic pressure for ε = 1/15 at x = L/4

Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .

 1011121314 Fig. 10. EHL pressure with transverse patterns at x 0 2 = 0, ε = 1/20

Fig. 15 .Fig. 16 .Fig. 17 .

 151617 Fig. 15. EHL average saturation with transverse patterns in the whole domain

Fig. 19 .

 19 Fig. 18. EHL pressure with longitudinal patterns at x 0 1 = -0.4

Fig. 20 .

 20 Fig. 20. Transverse roughness effects over the pressure in purely hydrodynamic and elastohydrodynamic cases at x 0 2 = 0

Fig. 21 .Fig. 22 .Fig. 23 .

 212223 Fig. 21. Transverse roughness effects over the saturation in purely hydrodynamic and elastohydrodynamic cases at x 0 2 = 0

Fig. 24 .

 24 Fig. 24. Influence of the piezoviscosity over the (average) deformation


Guy Bayada

Sébastien Martin

Carlos Vázquez

About a generalized Buckley-Leverett equation and lubrication multifluid flow

In this paper, we analyse the asymptotic system corresponding to a thin film flow with two different (non-miscible) fluids, from theoretical and numerical point of view. We also compare this model to the Elrod-Adams one, which is the reference model in tribology, as cavitation phenomena occur.

. Let be Ω =]0, L[ and let us denote by ∂Ω = {0, L} its boundary, by Q T the set ]0, T [×Ω and by Σ T the set ]0, T [×∂Ω. We introduce the ratio ε = µ g /µ l , µ l (resp. µ g ) being the viscosity of a fluid in liquid (resp. gaseous) phase. In view of cavitation-related phenomena, the fluid is

Introduction

The asymptotic behaviour of a single incompressible flow between two close surfaces in relative motion is described by the well-known Reynolds equation

∇ • h 3 6µ l ∇p = v 0 ∂ ∂x 1 h ,
in which h is the small gap between the two surfaces, v 0 the relative velocity of the surfaces, µ l the fluid viscosity and p the fluid pressure. This equation can be rigorously deduced from (Navier) Stokes system by means of an asymptotic analysis [START_REF] Bayada | The transition between the Stokes equations and the Reynolds equation: a mathematical proof[END_REF]. However, in some applications, the lubricant cannot be considered as a single fluid and a multifluid approach has to be introduced. For example, this happens when one of the surfaces has to be particularly protected from contact from the other one and it is covered by a specific fluid; this can be also modelled by the existence of a surface layer with a viscosity which is different from the one of the bulk fluid. Another phenomenon which falls into the scope of the multifluid approach is linked to cavitation, which introduces the existence of air bubbles inside the bulk fluid. Nevertheless, most of the multifluid problems in the lubrication area assume that the boundary between the two immiscible fluids is known [START_REF] Bayada | About thin film micropolar asymptotic equations[END_REF][START_REF] Jean Paulin | About the derivation of Reynolds law from Navier-Stokes equation for two non-miscible fluids[END_REF][START_REF] Tichy | Non-Newtonian lubrication with the convective Maxwell model[END_REF]. This assumption allows to obtain a slightly modified Reynolds equation. However, the real problem is a free boundary one and the position of the interface is defined by an additional unknown function. Using a three dimensional multifluid approach introduced by Nouri, Poupaud and Demay [START_REF] Nouri | An existence theorem for the multi-fluid Stokes problem[END_REF], a limit system has been derived by Bayada, Sabil and Paoli [START_REF] Sabil | Etude théorique et numérique d'un écoulement diphasique non-miscible et incompressible dans un domaine mince[END_REF][START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: from Stokes equations to Buckley-Leverett equation and Reynolds law[END_REF] with an asymptotic approach. This system describes the behaviour of the pressure and the relative saturation of the two fluids contained in a thin domain. However, this derivation is based upon an assumption on the shape of the interface. Thus, a mathematical study of the obtained asymptotic system has to be made in order to ensure its well-posedness. So far, this kind of result has been given only for the case when the surfaces surrounding the fluid are fixed (v 0 = 0), corresponding for example to the injection of a fluid through a fixed gap. Whether the value of v 0 is zero or not, the system consists of two equations : a generalized Buckley-Leverett equation and a generalized Reynolds lubrication equation.

However, the assumption of zero value for v 0 is not realistic for most of the lubrication problems in which the fluid is sheared, due to the difference of velocities between the surrounding surfaces. Thus, it is the purpose of this chapter to give an existence and uniqueness result for the asymptotic system with non zero value for the shear velocity. The main difficulty comes from the study of the generalized Buckley-Leverett equation. More precisely, taking the fact that v 0 is different from 0 prevents us from using the classical results about first order hyperbolic equations: in fact, the flux function is not autonomous and we have to guarantee that the saturation lies in the interval [0, 1], although the maximum principle is not a priori guaranteed anymore. This chapter is organized as follows:

Section 2 deals with the governing equations of the asymptotic system, obtained from the multifluid Stokes system. Thus, the generalized Buckley-Leverett / Reynolds system is presented, along with the physical assumptions related to realistic modellings. Section 3 is devoted to the analysis of the generalized Buckley-Leverett equation. Thus, we present the definition of a weak entropy solution of a scalar conservation law on a bounded domain, and we give some stability results (in particular, the saturation is a function with values in [0, 1]), along with an existence and uniqueness theorem, by using the concept of "semi Kružkov entropy-flux pairs". Moreover, we use a numerical scheme that allows to approximate the unique (physical) solution as the mesh size tends to 0. In Section 4, we present the analysis of the generalized Reynolds equation. In particular, we state an existence and uniqueness result for the pressure, and also a priori estimates in the H 1 or L ∞ norm which do not depend on the ratio ε of the viscosities. Section 5 deals with numerical computations. In particular, we present some numerical tests showing the importance of the shear effects on the saturation, pressure and velocity profiles but also on the boundary conditions. Then, we focus on cavitation phenomena : we show how our approach allows us to give some comprehensive details on the way the bifluid behaves in thin films. In particular, we compare the solution of the generalized Buckley-Leverett / Reynolds model to the solution of the Elrod-Adams model, which is a frequently used model in mechanical studies. This states that the Elrod-Adams model (which is heuristic) may be justified, at least numerically, by the bifluid model with an appropriate set of data.

supposed to be a lubricant: thus, µ l ≡ µ, the liquid phase lubricant being considered as a reference fluid, and the gaseous phase lubricant may be considered as air or gas. In that configuration, typically ε ∼ 10 -3 . Now we introduce the main equations: Generalized Buckley-Leverett equation: the saturation s is governed by a scalar conservation law:

∂ ∂t (h(x) s(t, x)) + ∂ ∂x (Q in (t) f (s) + v 0 h(x) g(s)) = 0, (t, x) ∈ Q T , (2.1)
where h is the normalized gap between the surfaces, Q in is the flow input to be related to the imposed saturation (corresponding to injection of the bifluid), v 0 is the shear velocity corresponding to the speed of the lower surface, and s denotes the reference (liquid) fluid saturation. The functions f and g, defining the flux, are described later. However, we point out the fact that f represents the classical contribution to the Buckley-Leverett flux, while g is a non-classical contribution induced by shear effects. Equation (2.1) is completed with the following initial and boundary conditions

s(0, •) = s 0 , on Ω (2.2) "s = s 1 ", on Σ T (2.3)
where the sense of the boundary condition (2.3) will be precised later. A common normalization of the flow input may be used:

Q in = θ in v 0 h(0)
, where θ in can be chosen equal to the boundary saturation s 1 . However, the flow input is mainly controlled by the input velocities of both fluids and the previous choice consists in neglecting the contribution of the fluid with the lowest viscosity (i.e. gas, by contrast with liquid).

This initial boundary value problem is weakly coupled with the following problem:

Generalized Reynolds equation: for a given saturation s, the pressure p obeys the following law:

∂ ∂x A(s) h 3 6µ l ∂p ∂x = v 0 ∂ ∂x (B(s)h) , (t, x) ∈ Q T , (2.4)
with the boundary condition:

p = 0, (t, x) ∈ Σ T . (2.5)
Let us mention that expressions of functions f , g, A and B are fully detailed in the Appendix (at the end of this paper). In [START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: from Stokes equations to Buckley-Leverett equation and Reynolds law[END_REF], the derivation of the coupled problem has been done under the assumption that the free boundary, which separates both phases, is a function belonging to L ∞ ((0, T ); BV (Ω)). Thus, this assumption (on the shape and on the regularity of the free boundary) prevents us from considering multi-layer flows, although they are supposed to be relevant (see, for instance, [START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: from Stokes equations to Buckley-Leverett equation and Reynolds law[END_REF]). The main reason for restricting ourselves to this particular type of free boundaries lies in the difficulty to compute an explicit expression of f and g otherwise. Therefore, we will restrict ourselves to cases (i) and (ii) (in which f , g, A and B can be fully computed):

(i) The reference (liquid) fluid is adhering to the lower (moving) surface (see Fig. 1).

(ii) The reference (liquid) fluid is adhering to the upper (fixed) surface (see Fig. 2).

We will discuss in Section 5 whether the choice of each assumption is relevant or not. Notice that f , g, A and B highly depend on the ratio ε. As it will be pointed out further, the shape of the flux functions and coefficients remains the same for both cases, as described by the forthcoming assumptions.

x z s(t, x)h(x) liquid gas h(x) L 0 We consider the following assumptions on the data: Assumption 1 (initial and boundary functions)

(i) s 0 ∈ L ∞ (Ω; [0, 1]), (ii) s 1 ∈ L ∞ (R + ; [0, 1]).
Assumption 2 (flow, shear velocity)

(i) Q in ∈ C 0 ([0, +∞]), (ii) ∃ Q min , Q max , 0 < Q min ≤ Q in ≤ Q max , (iii) v 0 > 0.
Assumption 3 (gap between the surfaces)

(i) h ∈ C 1 (R), (ii) ∃ h min , h max , 0 < h min ≤ h ≤ h max . Assumption 4 (auxiliary flux functions) (i) f ∈ C 1 ([0, 1]), f (0) = 0, f (1) = 1, f is a non-decreasing function, (ii) g ∈ C 1 ([0, 1]), g(0) = g(1) = 0,
Our purpose is to state an existence and uniqueness result for problem (2.1)-(2.5). In fact, the main difficulty is to state an existence and uniqueness result for the generalized Buckley-Leverett equation (2.1)-(2.3). Moreover, since s denotes a saturation, we have also to state that the (possible) solution takes its values in [0, 1]. Indeed, let us recall that the derivation of the generalized Buckley-Leverett equation is not fully rigorous in the sense that a strong assumption on the free boundary shape has been used. Let us recall also that the study of the generalized Buckley-Leverett equation including the shear term has been omitted by Paoli in [START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: from Stokes equations to Buckley-Leverett equation and Reynolds law[END_REF]. More precisely, we will prove the following result: Theorem 2.1 Under Assumptions 1-4, the bifluid model (2.1)-(2.5) admits a unique solution (p, s) ∈ L ∞ ((0, T ;

H 1 0 (Ω))× L ∞ (Q T ; [0, 1]
) in a sense that will be discussed later.

This will be a consequence of Theorems 3.9 and 4.1.

The generalized Buckley-Leverett equation

In a first subsection, we introduce an auxiliary problem and a corresponding "weak entropy solution", whose framework lies in the theory of scalar conservation laws on bounded domains. After establishing some results (existence, uniqueness, stability) on the properties of the auxiliary problem, we will prove, in a second subsection, how it is possible to reduce the generalized Buckley-Leverett problem to the auxiliary one. Finally, in the third subsection, we propose a numerical scheme whose solution converges to the "physical" solution (i.e. the weak entropy solution).

An auxiliary problem

Let us consider the following assumption:

Assumption 5
The "auxiliary boundary / initial data" satisfy :

(i) u 0 ∈ L ∞ ((0, 1); [0, 1]), (ii) u D ∈ L ∞ (R + ; [0, 1]).
The "auxiliary gap" satisfies:

(i) H ∈ C 1 (R 2 ), (ii) ∃ H min , H max , 0 < H min ≤ H ≤ H max .
We introduce the following scalar conservation law:

∂u ∂t + ∂ ∂x (f (u) + H(t, x) g(u)) = 0, (t, x) ∈]0, T [×]0, 1[ (3.1)
where f , g have been already defined and T > 0. Equation (3.1) is completed with the following initial and boundary conditions

u(0, •) = u 0 , on ]0, 1[, (3.2)
"u = u D ", on ]0, T [×∂]0, 1[. (3.3)
Existence and uniqueness of a solution for scalar conservation laws on unbounded domains has been solved in the pioneering work of Kružkov [START_REF] Kružkov | First order quasilinear equations with several independent variables[END_REF] who introduced the concept of weak entropy solution and related "Kružkov entropy-flux pairs". When dealing with bounded domains, introducing boundary conditions must be understood in a particular way: in the bounded domain setting, Bardos, Le Roux and Nédélec [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF] also proved existence and uniqueness of a weak entropy solution satisfying a "Kružkov entropy-flux pair" formulation which includes boundary terms, under some regularity assumptions on the data. In particular, the way the boundary condition is satisfied is known as the "BLN" condition. Nevertheless, this formulation would not be sufficient by two reasons:

1) boundary and initial conditions lack regularity in comparison to the framework of [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF], 2) it does not provide a stability result (we recall that we have to state that the possible solution is a function with values in [0, 1]).

The notion of weak entropy solution in the L ∞ framework is essentially due to Otto [START_REF] Otto | Initial-boundary value problem for a scalar conservation law[END_REF] who introduced the so-called "boundary entropy-flux pairs". A more complete exposition appears in [START_REF] Málek | Weak and measure-valued solutions to evolutionary PDEs[END_REF]. In fact, as it will be pointed out further, it is equivalent to use the "semi Kružkov entropy-flux pairs" and the "boundary entropy-flux pairs", at least in the case of a scalar conservation law with an autonomous flux. In the case of scalar conservation laws with non-autonomous fluxes, a "boundary entropy-flux pairs" formulation is not obvious anymore, but using the concept of "semi Kružkov entropy-flux pairs" [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF][START_REF] Serre | Systèmes de lois de conservation[END_REF][START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF][START_REF] Ammar | Scalar conservation laws with general boundary condition and continuous flux function[END_REF]

1 0 e T 0 (u -k) ± ∂ϕ ∂t + Φ ± [f ] (u, k) + H Φ ± [g] (u, k) ∂ϕ ∂x -sgn ± (u -k) ∂H ∂x g(k) ϕ dx dt + 1 0 (u 0 (x) -k) ± ϕ(0, x) dx + M e T 0 (u D (t, 1) -k) ± ϕ(t, 1) dt + M e T 0 (u D (t, 0) -k) ± ϕ(t, 0) dt ≥ 0, (3.4
)

for all k ∈ [0, 1], ϕ ∈ D([0, T [×[0, 1]
) and ϕ ≥ 0, the constant M being defined by

M = f ′ L ∞ (0,1) + k max g ′ L ∞ (0,1) . (3.5)
The functions u → (u -k) ± are the so-called "semi Kružkov entropies" (see [START_REF] Carrillo | Entropy solutions for nonlinear degenerate problems[END_REF][START_REF] Serre | Systèmes de lois de conservation[END_REF][START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF]), defined by

(u -k) + = u -k, if u ≥ k, 0, otherwise. and (u -k) -= (k -u) + .
The functions Φ ± [f ] and Φ ± [g] are the corresponding "semi Kružkov fluxes" defined by

Φ ± [f ] (u, k) = sgn ± (u -k)(f (u) -f (k)), Φ ± [g] (u, k) = sgn ± (u -k)(g(u) -g(k)),
where u → sgn ± (u) is the derivative of the function u → u ± with value 0 at point 0.

The functions Φ ± [f ] (•, k) and Φ ± [g] (•, k) are defined on [0, 1] (as f and g). Thus, both definitions of solutions (weak entropy and strong solutions) make sense only if the function u takes values in [0, 1] a.e., which is not supposed a priori in Definition 1. However, the set is preserved as shown by the following proposition. [f ] , Φ ± [g] by the functions f and g, Φ ±

[f ] and Φ ± [g] defined by:

f (v) =    0, if v < 0 f (v), if 0 ≤ v ≤ 1 1, if v > 1 , Φ ± [f ] (v, k) = sgn ± (v -k)( f (v) -f (k)), g(v) =    0, if v < 0 g(v), if 0 ≤ v ≤ 1 0, if v > 1 , Φ ± [g] (v, k) = sgn ± (v -k)( g(v) -g(k)),
and if u is a weak entropy solution of problem (3.1)-(3.3), then 0 ≤ u ≤ 1 a.e. on ]0, T [×]0, 1[.
Proof First let us notice that the following properties hold:

Φ ± [f ] (v, k) ≤ f ′ L ∞ (0,1) (v -k) ± , Φ ± [g] (v, k) ≤ g ′ L ∞ (0,1) (v -k) ± .
Now, set k = 0 in Inequality (3.4). Since we have (u 0 ) -= 0, (u D ) -= 0 (see Assumption 5), g(0) = 0 (see Assumption 4), the three last terms in (3.4) vanish. Thus, we have

1 0 e T 0 u -∂ϕ ∂t + Φ - [f ] (u, 0) + H Φ - [g] (u, 0) ∂ϕ ∂x dx dt ≥ 0. (3.6) Now, let (τ, R) be such that 0 < τ < T , δ = T -τ and 0 < R < 1. Let r ∈ D(R + ) be such that: r is non-increasing, r ≡ 1 on [0, R + Mτ ], r ≡ 0 on [R + Mτ + δ/2, +∞). Then, choosing ϕ(t, x) = τ -t τ χ (0,τ ) (t) r(x + Mt) in Inequality (3.6) leads to - 1 τ 1 0 τ 0 u -r(x + Mt) dt dx + 1 0 τ 0 τ -t τ r ′ (x + Mt) M u -+ Φ - [f ] (u, 0) + H Φ - [g] (u, 0) ≥ 0. (3.7) Now, since Φ - [f ] (u, 0) ≤ f ′ L ∞ (0,1) u -, Φ - [g] (u, 0) ≤ g ′ L ∞ (0,1)
u -, we use Equality (3.5) to obtain:

M u -+ Φ - [f ] (u, 0) + H Φ - [g] (u, 0) ≥ 0. Moreover, since r ′ (x + Mt) ≤ 0, the second term of the left-hand side of (3.7) is negative. Now, since r(x + Mt) = 1, for all (t, x) ∈ (0, τ ) × (0, R) and r ≥ 0, the left-hand side of (3.7) is upper bounded by - 1 τ R 0 τ 0 u -dx dt
which is consequently non-negative. Therefore, we have u -= 0 on (0, τ ) × (0, R). Now, passing to the limit with R → 1 and τ → T , we have u ≥ 0 a.e. Similarly, by choosing k = 1 in Inequality (3.4) (with the semi entropies u → (u -1) + ), we prove that u ≤ 1 a.e. on ]0, T [×]0, 1[. Remark 3.2 Interestingly, the fact that the flux function in the auxiliary problem (3.1)-(3.3) is not autonomous involves a major difference with the autonomous case, concerning the stability intervals. Indeed, the following properties can be easily shown:

- This is due to the properties of g, in particular g(0) = g(1) = 0.

if
Remark 3.3 A function u which satisfies Definition 1 is a weak solution in a classical sense. Indeed, for every ϕ ∈ H 1 0 (]0, T [×]0, 1[), we write ϕ = ϕ + -ϕ -, with ϕ + = max(ϕ, 0) and ϕ -= -min(ϕ, 0); obviously, ϕ ± ∈ H 1 0 (]0, T [×]0, 1[);
thus adding the two inequalities (corresponding to each "semi Kružkov entropy-flux pair") gives:

e T 0 1 0 u -k ∂ϕ ± ∂t + sgn(u -k) (f (u) -f (k)) + H (g(u) -g(k)) ∂ϕ± ∂x - ∂H ∂x sgn(u -k) g(k) ϕ ± dx dt ≥ 0. Now, taking k = ± u L ∞ (]0, e T [×]0,1[) gives: e T 0 1 0 u ∂ϕ ± ∂t + (f (u) + H g(u)) ∂ϕ ± ∂x dx dt = 0,
Now, by means of substraction, we immediatly obtain Let us explain the way the boundary condition is satisfied:

Definition 2 (Boundary entropy-flux pairs, [START_REF] Otto | Initial-boundary value problem for a scalar conservation law[END_REF]) Let us define the flux f :

(t, x, s) -→ f (s) + H(t, x) g(s). The pair ( Ĥ, Q) belonging to C 1 (R 2 ) × C 1 (]0, T [×]0, 1[×R 2
) is said to be a "boundary entropy-flux pair" (for the flux f ) if it satisfies:

1. for all w ∈ R, s → Ĥ(s, w) is a convex function, 2. ∀(t, x) ∈]0, T [×]0, 1[, ∀w ∈ R, ∂ 3 Q(t, x, s, w) = ∂ 1 Ĥ(s, w)∂ 3 f (t, x, s), 3. ∀w ∈ R, Ĥ(w, w) = 0, Q(•, •, w, w) = 0, ∂ 1 Ĥ(w, w) = 0. Proposition 3.4 (Boundary condition, [16]) Let u ∈ L ∞ (]0, T [×]0, 1[) be a weak entropy solution of problem (3.1)-(3.3). Then, ess lim ̺→0 + e T 0 Q(t, 1, u(t, 1 -̺), u D (1)) β(t, 1) -Q(t, 0, u(t, ̺), u D (0)) β(t, 0) dt ≥ 0,
for all "boundary entropy-flux pairs" ( Ĥ, Q), ∀β ∈ L 1 (Σ T ), β ≥ 0 a.e. (3.8)

Remark 3.5 Now let us give some comprehensive details on the way to understand the boundary condition. This has been given in [START_REF] Málek | Weak and measure-valued solutions to evolutionary PDEs[END_REF][START_REF] Otto | Initial-boundary value problem for a scalar conservation law[END_REF][START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF]: in general, the problem should be overdetermined and the boundary equality cannot be required to be assumed at each point of the boundary, even if the solution is a regular function. But, with additional assumptions, the more comprehensive "BLN" condition is recovered:

(i) If u admits a trace, i.e. there exists u |b ∈ L ∞ (]0, T [×∂]0, 1[) such that ess lim ̺→0 + e T 0 |u(t, 1 -̺) -u |b (t, 1)| + |u(t, ̺) -u |b (t, 0)| dt = 0,
then, Inequality (3.8) is equivalent to the following inequality (see [START_REF] Dubois | Boundary conditions for nonlinear hyperbolic systems of conservation laws[END_REF][START_REF] Otto | Initial-boundary value problem for a scalar conservation law[END_REF])

Q(•, 1, u |b (•, 1), u D (•, 1)) ≥ 0, a.e. on ]0, T [, (3.9)
Q(•, 0, u |b (•, 0), u D (•, 0)) ≤ 0, a.e. on ]0, T [. (3.10)
Moreover, considering the particular "boundary fluxes"

H + δ (z, w) = (max(z -w, 0)) 2 + δ 2 1/2 -δ, (3.11)
Q + δ (t, x, z, w) = z w ∂ 1 H + δ (λ, k) ∂ 3 f (t, x, λ) dλ, (3.12)
and

H - δ (z, w) = (min(w -z, 0)) 2 + δ 2 1/2 -δ, (3.13)
Q - δ (t, x, z, w) = z w ∂ 1 H - δ (λ, k) ∂ 3 f (t, x, λ) dλ, (3.14)
and letting δ → 0, we have the following uniform convergences:

Q ± δ (t, x, s, w) → sgn ± (s -w) f (t, x, s) -f (t, x, w)) .
Finally, taking the following "boundary flux" in Inequalities (3.9) and (3.10)

Q(t, x, s, w) = sgn + (s -max(w, k)) f (t, x, s) -f (t, x, max(w, k))) +sgn -(s -min(w, k)) f (t, x, s) -f (t, x, min(w, k)))
yields to the classical BLN condition given by Bardos, Le Roux and Nédélec [START_REF] Bardos | First order quasilinear equations with boundary conditions[END_REF], that is:

sgn(u |b (t, 1) -u D (t, 1))( f (t, 1, u |b (t, 1)) -f (t, 1, k)) ≥ 0, (3.15) sgn(u |b (t, 0) -u D (t, 0))( f (t, 0, u |b (t, 0)) -f (t, 0, k)) ≤ 0, (3.16) for a.e. (t, r) ∈]0, T [×∂]0, 1[, ∀k ∈ [min(u |b , u D ), max(u |b , u D )].
(ii) Assume that u admits a trace on the boundary, then Inequalities (3.15) and (3.16) can be simplified in the following cases:

• If ∂ 3 f is a positive function on ]0, 1[, then u(•, 0) = u D (and nothing is imposed at x = 1, i.e. the condition u(•, 1) = u D is not active). • If ∂ 3 f is a negative function on ]0, 1[, then u(•, 1) = u D (and nothing is imposed at x = 0, i.e. the condition u(•, 0) = u D is not active).
Thus the boundary conditions may be "active" only on a part of the boundary. Unfortunately, in the case of the flux f with Assumptions 4 and 5, monotonicity with respect to the third variable lacks, and we have to deal with "relaxed" Dirichlet boundary conditions. Now, we conclude this subsection with the main result related to the auxiliary problem: Proof Existence is proved using the (classical) parabolic approximation, which consists of adding an artificial diffusive term in the right-hand side of the hyperbolic equation (vanishing viscosity method). Next, passing to the limit on the diffusive parameter gives the existence result. Uniqueness is obtained from the Kružkov method of doubling variables. The complete proof is, in fact, valid for first order quasilinear equations in any space dimension, and completely detailed in [START_REF] Martin | First order quasilinear equations with boundary conditions in the L ∞ framework[END_REF]. Let us sketch the main steps of the proof:

Existence of the solution:

• Step 1 -The parabolic approximation:

∂u ε ∂t + ∂ ∂x f (u ε ) + H(t, x) g(u ε ) = ε ∂ 2 u ε ∂x 2 , on ]0, T [×]0, 1[, (3.17)
u ε (0, •) = u 0 ε , on ]0, 1[, (3.18) u = u D ε , on ]0, T [×∂]0, 1[, (3.19)
where

u D ε ∈ C 2 (]0, T [×∂]0, 1[; [0, 1]) and u 0 ε ∈ C 2 ([0, T ] × [0, 1]; [0, 1]
) satisfy compatibility conditions on [0, T ] × ∂]0, 1[. Notice that the initial / boundary conditions need to be regularized.

• Step 2 -Entropy formulation of the parabolic problem: Introducing a regularizing function 1 , let (u, u D , u 0 ) satisfy equations (3.17)-(3.19), with regular boundary / initial conditions (subscripts are dropped for convenience). Then, for all ϕ ∈ D(] -∞, T [×R), 1 The following definition is introduced:

Definition 3 Let µ be a sufficient small positive number, and let us define the function:

s(x) =  min (dist(x, ∂]0, 1[), µ) , if x ∈]0, 1[, -min (dist(x, ∂]0, 1[), µ) , if x ∈ R\]0, 1[. Let ξε be defined by ξε(x) = 1 -exp " - M + εR ε s(x) « , with R = sup 0<s(x)<µ |∆s(x)|.
for all k ∈ R,

e T 0 1 0 (u -k) ± ∂ϕ ∂t + sgn ± (u -k) (f (u) -f (k)) + H (g(u) -g(k)) ∂ϕ ∂x -sgn ± (u -k) ∂H ∂x f (k) ϕ + ε (u -k) ± ∂ 2 ϕ ∂x 2 ξ ε dx dt (3.20) + Ω (u 0 -k) ± ϕ(0, x) ξ ε dx ≥ -2ε e T 0 1 0 (u -k) ± ∂ 2 ϕ ∂x 2 ∂ 2 ξ ε ∂x 2 dx dt -(M + Rε) e T 0 (u D (t, 1) -k) ± ϕ(t, 1) dt + e T 0 (u D (t, 0) -k) ± ϕ(t, 0) dt , • Step 3 -L 1 -stability result: Let (u 1 , u D 1 , u 0 1 ), (u 2 , u D 2 , u 0
2 ), satisfy equations (3.17)-(3.19), the boundary / initial conditions being regular. Then, for all t ∈ (0, T ),

1 0 u 1 (t, •) -u 2 (t, •) ξ ε dx ≤ 1 0 u 0 1 -u 0 2 ξ ε dx + (M + Rε) e T 0 u D 1 (•, 1) -u D 2 (•, 1) + u D 1 (•, 0) -u D 2 (•, 0) dt .
• Step 4 -BV estimates: Let (u, u D , u 0 ) satisfy equations (3.17)- (3.19), the corresponding boundary / initial conditions being regular. We suppose furthermore that u D has a smooth extension to Q T , denoted u D . Then, there exists a constant λ which only depends on u 0 Ω , u D ΣT , T , Ω =]0, 1[, f , g and H such that

sup t∈(0, e T ) Ω ∂u ∂t (t, •) + ∂u ∂x (t, •) ≤ λ, (3.21)
where we have used the notation

u 0 Ω = 1 0 ∂ 2 u 0 ∂x 2 + ∂u 0 ∂x + u 0 dx, u D ΣT = sup ]0, e T [×]0,1[ ∂ 2 u D ∂x 2 + ∂u D ∂t + ∂u D ∂x + u D + e T 0 1 0 ∂ 2 ∂x 2 ∂u D ∂t + ∂ 3 u D ∂x 3 + ∂ 2 u D ∂t 2 + ∂ ∂x ∂u D ∂t + ∂ 2 u D ∂x 2 .
• Convergence of the sequence: Let u ε be the unique solution of Equations (3.17)- (3.19) corresponding to (regular) initial / boundary conditions (u 0 ε , u D ε ) and let lim

ε→0 u D ε = u D in L 1 (]0, T [×∂]0, 1[), lim ε→0 u 0 ε = u 0 in L 1 (]0, 1[),
where

u D ∈ L ∞ (]0, T [×∂]0, 1[; [0, 1]) and u 0 ∈ L ∞ (]0, 1[; [0, 1]). Then, the sequence {u ε } ε converges to some function u ∈ L ∞ (]0, T [×]0, 1[; [0, 1]) in C 0 ([0, T ], L 1 (]0, 1[)).
Moreover u is a weak entropy solution of problem (3.1)- (3.3). This result is obtained as follows: our goal is to let ε tend to 0 in Equations (3.17)- (3.19). Nevertheless, we cannot apply the BV estimates on u ε because u D ε , u 0 ε satisfy compatibility conditions but do not necessarily have an extension over [0, T ] × [0, 1] with sufficient regularity. Thus, we introduce, by means of construction, (u D ε,h , u 0 ε,h ) which both satisfy compatibility conditions and have an extension over [0, T ] × [0, 1] with sufficient regularity. Moreover, (u D ε,h , u 0 ε,h ) are uniformly "close" to (u D ε , u 0 ε ) (as h → 0, uniformly with respect to ε), which implies that u ε,h is "close" to u ε (due to the L 1 -stability result). Then, we apply Arzelà-Ascoli theorem to the sequence {u ε } in order to prove that it is relatively compact in C 0 ([0, T ]; L 1 (]0, 1[)). Of course, we have to verify that the sequence satisfies the hypotheses of the theorem (equicontinuity and pointwise relative compactness): for this, we use the properties of u ε,h and the fact that u ε,h is "close" to u ε,h . The limit function is an entropy solution of problem (3.17)-(3.19) by passing to the limit in the entropy formulation of the parabolic problem.

Uniqueness:

• Step 1 -Comparison: Let u ∈ L ∞ (]0, T [×]0, 1[) (resp. v ∈ L ∞ (]0, T [×]0, 1[)) be an entropy solution of problem (3.1)-(3.
3) with initial and boundary conditions (u

0 , u D ) ∈ L ∞ (]0, 1[) × L ∞ (]0, T [×∂]0, 1[) (resp. (v 0 , v D ) ∈ L ∞ (]0, 1[) × L ∞ (]0, T [×∂]0, 1[)); then - e T 0 1 0 u -v ∂β ∂t + sgn(u -v) (f (u) -f (v)) + H (g(u) -g(v)) ∂β ∂x ≤ M e T 0 u D (t, 1) -v D (t, 1) β(t, 1) + u D (t, 0) -v D (t, 0) β(t, 0) dt + 1 0 u 0 (x) -v 0 (x) β(0, x) dx for all β ∈ D((-∞, T )×R d ).
The result is obtained using Kružkov's method of doubling variables, adapted to the entropy formulation, and a specific treatment of the boundary terms. • Step 2 -Uniqueness: Considering the previous integral inequality with v D = u D and v 0 = u 0 and a test function which only depends on time t, we get:

e T 0 1 0 u -v α ′ (t) dx dt ≥ 0, (3.22)
for all α ∈ D(-∞, T ). Then, for an interval [t 0 , t 1 ] ⊂]0, T [, we can use in Inequality (3.22) the characteristic function of [t 0 , t 1 ], properly mollified, and pass to the limit on the mollifier parameter:

Ω u(t 1 , •) -v(t 1 , •) ≤ Ω u(t 0 , •) -v(t 0 , •)
As t 0 tends to 0, and using the fact that v 0 = u 0 , uniqueness is straightforward.

Existence and uniqueness of a weak entropy solution for the Buckley-Leverett problem

In this subsection, we show that problem (2.1)-(2.3) can be reduced to an auxiliary problem as the one described in the previous subsection, namely problem (3.1)-(3.3). Then, it suffices to use the results established for the auxiliary problem.

Definition 4 (Direct reduction) Let us consider the following changes of variables:

Y (x) = x 0 h(t) dt L 0 h(x) dx , T (t) = L t 0 Q in (s) ds L 0 h(x) dx .
We also define the inverse functions of Y and T , respectively denoted Y -1 and T -1 .

Definition 5 (Inverse reduction) Let Y -1 and T -1 be defined as the respective unique solution of the following Cauchy problems: Proof First, we define

   dY -1 dy (y) = 1 h(Y -1 (y)) L 0 h(x) dx, Y -1 (0) = 0,    dT -1 dτ (τ ) = 1 L Q in (T -1 (τ )) L 0 h(x) dx, T -1 (0) = 0.
u(T (t), Y (x)) = s(t, x), (3.23)
so that Equation (2.1) reduces to

∂ ∂τ u(τ, y) + ∂ ∂y (f (u) + H(τ, y) g(u)) = 0, (τ, y) ∈ (0, T ) × (0, 1) (3.24) with H(τ, y) = v 0 (Q in • T -1 ) (τ ) h • Y -1 (y).

About a generalized Buckley-Leverett equation and lubrication multifluid flow 15

The initial and boundary conditions are modified as follows:

u(0, y) = u D (0, y) = s 0 • Y -1 (y), y ∈]0, 1[ (3.25) "u(τ, y) = u D (τ, y) = s 1 T -1 (τ ), Y -1 (y) ", (τ, y) ∈ (0, T ) × ∂]0, 1[ (3.26)
As a consequence, turning back to the original variables immediatly gives:

Definition 6 A function s ∈ L ∞ (Q T ; [0, 1]
) is said to be a weak entropy solution of problem (2.1)-( 2.3) if it satisfies

QT h(x) (s -k) ± ∂ϕ ∂t + Q in (t) Φ ± [f ] (s, k) + v 0 h(x) Φ ± [g] (s, k) ∂ϕ ∂x -v 0 sgn ± (s -k) h ′ (x) g(k) ϕ dx dt + Ω h(x) (s 0 (x) -k) ± ϕ(0, x) dx + L T 0 (s 1 (t, 1) -k) ± ϕ(t, 1) dt + L T 0 (s 1 (t, 0) -k) ± ϕ(t, 0) dt ≥ 0 (3.27) for all k ∈ [0, 1], ϕ ∈ D([0, T [×R), ϕ ≥ 0.
Here, we can choose:

L = Q max max( f ′ L ∞ (0,1) ) + v 0 h max Q max Q min max( g ′ L ∞ (0,1) ).
Theorem 3.9 Under Assumptions 1-4, problem (2.1)-(2.3) admits a unique weak entropy solution.

Numerical scheme for the generalized Buckley-Leverett equation

We study the numerical method simulating the solution of problem (3.24)-(3.26) (for convenience) which is equivalent to problem (2.1)-(2.3) up to a change of variables.

The following theorem gives a convergence result which enables us to use some classical finite volume methods (such as the Lax-Friedrichs scheme, for instance) in order to compute the weak entropy solution.

Theorem 3.10 (Vovelle [START_REF] Vovelle | Convergence of finite volume monotone schemes for scalar conservation laws on bounded domains[END_REF]) Let us consider a finite volume scheme with monotone fluxes associated to problem (3.24)-(3.26) and its corresponding numerical solution u T ,k . Then (u T ,k ) strongly converges to the weak entropy solution u in L p loc (R + , Ω) for every p ∈ [2, +∞[. From a theoretical point of view, it appears that the simulation of the Buckley-Leverett problem can be easily done thanks to the earlier convergence result. Nevertheless, let us recall that f and g highly depend on the ratio of the viscosities, namely ε (see Fig. 14, 15 and 16 in the appendix of the paper). In practical situations, ε may be small (10 -3 for an air-liquid mixture). Thus, let us study the behaviour of the scheme for small values of ε. For this, let us consider a uniform admissible mesh, whose step size is denoted ∆x, with a time step ∆t. We first notice that:

• f ′ L ∞ and g ′ L ∞ tend to explode as ε tends to 0. More precisely, we have (after omitted computations, using Maple, based on an asymptotic expansion)

f ′ L ∞ = O ε -1/3 , g ′ L ∞ = O ε -1/3
, so that we get also M = O ε -1/3 (see Equation (3.5)). • Moreover, Id denoting the identity application, we have (using Maple again)

f L 1 = O ε 1/3 , Id -g L 1 = O ε 1/3 .
Then, ∆x needs to be adapted to the value of ε in order to describe phenomena related to the boundary layer, so that the mesh size should satisfy

∆x = O(ε 1/3 ) (3.28)
in order to describe boundary layers of f and g. Now, we recall that the following CFL condition has to be imposed in order to ensure the stability of the numerical scheme:

∃ξ ∈]0, 1[, ∆t ≤ (1 -ξ) ∆x M .
Thus, we obtain the order of the time step:

∆t = O(ε 2/3 ). (3.29)
Now, it clearly appears that it becomes difficult to simulate the Buckley-Leverett equation for too small values of ε. Indeed, for each time step, the number of computations increases with 1/∆x; similarly, the number of times steps increases with 1/∆t so that it becomes more and more difficult to attain a (possible) stationary solution by passing to the limit in the evolutive problem.

4 The generalized Reynolds equation

Existence and uniqueness

Definition 7 Let s ∈ L ∞ (Q T ; [0, 1]). A function p ∈ L ∞ ((0, T ); H 1 0 (Ω)) is a weak solu- tion of problem (2.4)-(2.5) if it satisfies L 0 A(s) h 3 6µ l ∂p ∂x ∂v ∂x dx = v 0 L 0 B(s)h ∂v ∂x dx, ∀v ∈ H 1 0 (Ω) (4.1)
for almost every t ∈ (0, T ).

About a generalized Buckley-Leverett equation and lubrication multifluid flow 17

Next, we establish the existence and uniqueness of solution for the generalized Reynolds equation:

Theorem 4.1 Under Assumptions 1-4, problem (2.4)-(2.5) admits a unique solution p in the sense of Definition 7. Moreover, we have the following estimates:

p L ∞ ((0,T );H 1 (Ω)) ≤ C 1 , p L ∞ (QT ) ≤ C 2 ,
uniformly with respect to ε.

Proof Obviously, time t plays the role of a parameter. Thus, considering a fixed time t, existence and uniqueness of the pressure follows from Lax-Milgram lemma applied to an elliptic problem. To achieve the proof, it is sufficient to state that the coefficients are bounded and that the A is coercive, which is obvious from the following properties

1 ≤ A(s) ≤ 1/ε, ∀s ∈ [0, 1], 0 ≤ B(s) ≤ 2, ∀s ∈ [0, 1],
which are easily obtained from the definitions of A and B (see the Appendix at the end of the paper). Here, ε denotes the ratio of the viscosities µ g /µ l . Thus, using p(t, •) (t ∈ (0, T )) as a test-function in the variational formulation (4.1), we have:

h 3 min L 0 ∂p ∂x (t, •) 2 ≤ L 0 A(s(t, •)) h 3 ∂p ∂x (t, •) 2 = 6v 0 µ l L 0 B(s(t, •)) h ∂p ∂x (t, •) ≤ 12v 0 µ l h max L 0 ∂p ∂x (t, •).
Next, using the Cauchy-Schwarz inequality, we have:

∂p ∂x (t, •) L 2 (Ω) ≤ 12v 0 µ l h max L 1/2 h 3 min .
The estimate in the H 1 norm is straightforward from Poincaré-Friedrichs inequality.

The estimate in the L ∞ norm comes from the fact that H 1 (Ω) ⊂ L ∞ (Ω) with compact injection. Since the earlier inequality holds for almost every t ∈ (0, T ), the proof is concluded.

Simulation of the generalized Reynolds equation

Suppose that the saturation s can be computed for each time step. Then, the pressure p is also obtained at each time step by using any numerical method related to an elliptic problem: finite difference discretization, shooting method with a Runge-Kutta solver (after some linear interpolation procedure on the saturation) or finite elements method.

5 Numerical simulation

Influence of the shear effects and boundary conditions

As pointed out in the previous sections, the shear effects play a crucial role on the analysis of the generalized Buckley-Leverett equation. Indeed, we have already mentioned that it leads to the non-autonomous property of the flux, and also to a relaxation of the Dirichlet conditions: without shear velocity, the boundary condition is active only at point x = 0 because of the monotonicity of f . But when the shear effects are included, this is no longer true due to the presence of the partial flux g. We are interested in the contribution of the shear term in the Buckley-Leverett equation. For this, we compare the following regimes:

• v

(1) 0 = 0: injection of a fluid through a fixed gap, • v

(2) 0 = 1: injection of a fluid through a fixed gap between two surfaces in relative motion (shear effects).

The main difference lies in the properties of the Buckley-Leverett flux function. For this, we choose the following numerical data:

• The lubricant is adhering to the moving surface (Case (i)),

• Geometrical data: Ω =]0, 1[, h(x) = (2x -1) 2 + 1 2 ,
• Reference viscosity:

µ l = µ g = 1,
• Boundary condition:

s 1 = θ in with θ in = 0.37, • Flow input: Q in = v (2)
0 θ in h(0), • The mesh grid has 600 elements,

• The CFL condition is given by ∆t ∆x M = 0.9.

We start from two different initial conditions:

s (1) 0 (x) = (θ in -0.01)(2x -1) 6 + 0.01, s (2)
0 (x) = (θ in -0.99)(2x -1) 6 + 0.99, (s

(i)
0 is called initial condition (i)) and we observe (see Fig. 3456) that in both regimes (v 0 = 0 or v 0 = 0), the numerical stationary solution (obtained for T = T f ) does not depend on the initial condition. Moreover, the shear effects involve a major difference with the autonomous case: the stationary solution is not constant but contains balanced effects due to the non-autonomous flux. Now, let us focus on the boundary conditions. As it was pointed out, the boundary conditions may be active only on some part of the boundary. Moreover, Theorem 3.10 provides a strong convergence result, up to a possible boundary layer at each point of the boundary. This kind of behaviour is illustrated on Fig. 7. As in the previous tests, the lubricant is adhering to the moving surface (Case (i)), the geometry is unchanged (converging-diverging profile) and the following numerical data have been considered:

• Shear velocity: v 0 = 1, • Viscosities: µ l = 1, µ g = 10 -3 µ l , • Flow input: Q in = v 0 h(0) g(θ in )/(1 -f (θ in )) with θ in = 0.385,
• The mesh grid has 1400 elements,

• The CFL condition is given by ∆t ∆x M = 0.9.

The initial condition is s 0 ≡ θ in and the boundary condition also takes the value s 1 = θ in (in the sense that has been precised before).

About a generalized Buckley-Leverett equation and lubrication multifluid flow 19

About cavitation phenomena in lubrication theory

Lubricated devices are generally made of two surfaces which are closely spaced, the annular gap being filled with some lubricant. The radial clearance is very small, typically ∆r/r = 10 -3 for infinite oil lubricated bearings wich allows us to use the Reynolds equation:

d dx h 3 6µ l dp dx = v 0 dh dx
where µ l is the lubricant viscosity, p the pressure distribution, h the height between the two surfaces and v 0 the relative speed of the surfaces.

Nevertheless, this modelling does not take into account cavitation phenomena: cavitation is defined as the rupture of the continuous film due to the formation of air bubbles and makes the Reynolds equation no longer valid in the cavitation area. In order to make it possible, various models have been used, the most popular perhaps being variational inequalities which have a strong mathematical basis but lack physical evidence. Thus, a more realistic model, the Elrod-Adams model, is often used, assuming that the cavitation region is a fluid-air mixture and introducing an additional unknown θ (the saturation of fluid in the mixture) (see [START_REF] Capriz | Partial lubrication of full cylindrical bearings[END_REF][START_REF] Coyne | Conditions for the rupture of a lubricating film, Part 1[END_REF][START_REF] Coyne | Conditions for the rupture of a lubricating film, Part 2[END_REF][START_REF] Elrod | A computer program for cavitation[END_REF]). The model includes a modified Reynolds equation, with the following formulation:

(P)    d dx h 3 6µ l dp dx = v 0 dθh dx , p ≥ 0, 0 ≤ θ ≤ 1, p(1 -θ) = 0,
with, for instance, the boundary conditions:

• Dirichlet conditions : p = 0, at x = L,

• Neumann conditions : v 0 θ h - h 3 6µ l dp dx = Q in , at x = 0.
Introducing the domains

Ω + = {x ∈ Ω, p(x) > 0} (lubricated region), Ω 0 = {x ∈ Ω, p(x) = 0} (cavitated region), Σ = ∂Ω + ∪ Ω (free boundary),
the free boundary Σ separates a full film area Ω + from a cavitated area Ω 0 . Notice that Q in denotes the flow input, which can be expressed as

Q in = v 0 θ in h(0), with 0 ≤ θ in ≤ 1.
This model has been studied in [START_REF] Bayada | Analysis of a free boundary problem in partial lubrication[END_REF]. In particular, Neumann conditions are proved to be compatible with homogeneous Dirichlet conditions if and only if the normalized flow input θ in belongs to some interval. More precisely, there exists (θ, θ) such that

• if 0 ≤ θ in ≤ θ, then the following obvious solution is obtained:

p := 0, θ := θ in h ,
• if θ ≤ θ in ≤ θ, then a "starvation" phenomenon occurs; the domain is divided into three areas: starvation area (Ω

0 in which p ≡ 0), then a satured area (Ω + in which p > 0) and then another cavitated area (Ω (2) 0 in which p ≡ 0).

• if θ in > θ, then the unique solution (p, θ) satisfies p(0) > 0 and p(L) = 0, which lacks physical evidence: in infinite journal bearings, both boundaries are located at the same place, corresponding to the supply groove and discontinuities of the pressure are not admissible. In fact, the Elrod-Adams is not valid for this type of regimes. Now, we want to compare the generalized Buckley-Leverett / Reynolds model to the Elrod-Adams model, which is motivated by the multifluid approach: indeed, the Elrod-Adams model describes partial lubrication of a device composed of a lubricant in two phases (liquid and gas) but strictly focuses on the liquid aspect. The generalized Buckley-Leverett / Reynolds model takes into account the influence of both phases. Thus, recalling that ε ∼ 10 -3 in the case of a bifluid composed of liquid/gas, we want to numerically compare the two models. More precisely, we aim at comparing the numerical stationary solution of the generalized Buckley-Leverett / Reynolds problem to the solution of the Elrod-Adams problem. For this, let us motivate the choice of the data:

• Geometrical data: Ω =]0, 1[, h(x) = (2x -1) 2 + 1 2 ,
• Shear velocity, reference viscosity (lubricant viscosity) : v 0 = 1 and µ l = 1,

• Flow input : θ in = 0.385. Let us notice that the same flow input Q in has been considered in the Buckley-Leverett / Reynolds model and the Elrod-Adams model.

The choice of such a flow input is designed to compare the Buckley-Leverett / Reynolds model to the Elrod-Adams model with the same data. With the chosen value, the Elrod-Adams solution contains starvation along with a well-known discontinuity in the saturation function corresponding to the rupture of the film thickness; thus, same phenomena observed with the Buckley-Leverett / Reynolds model would lead to some kind of a justification of the Elrod-Adams model. Unfortunately, it is impossible to obtain a Buckley-Leverett saturation with strictly value 1 at one point of the domain with such a choice, as it is stated in the following result:

Proposition 5.1 Let s be a stationary solution of the generalized Buckley-Leverett equations (2.1)-(2.3), and let us consider the following sentences:

(i) There exists x ∈ Ω such that s(x) = 1.

(ii) The relative flow input satisfies the so-called "full saturation condition":

θ in = g(θ in ) 1 -f (θ in ) .
Then we have (i) ⇒ (ii).

Proof Let us assume that s is a stationary solution of Equation (2.1). Then, after integration, Q in f (s(x)) + v 0 h(x) g(s(x)) = C, ∀x ∈ Ω, C being a constant. This constant is determined by the value of the earlier function at one boundary of Ω (notice that h(0) = h(1)), namely C = Q in f (θ in ) + v 0 h(0) g(θ in ). Since f (1) = 1 and g(1) = 0, if About a generalized Buckley-Leverett equation and lubrication multifluid flow 21 there exists x 0 ∈ Ω, such that s(x 0 ) = 1, one gets

Q in f (s(x 0 )) + v 0 h(x 0 ) g(s(x 0 )) = Q in . (5.1)
Thus, from the previous equations,

Q in = v 0 h(0) θ in = v 0 h(0) g(θ in ) 1 -f (θ in )
.

Proposition 5.2 We have the following properties:

• For all θ in ∈]0, 1[, if s is a stationary solution of the generalized Buckley-Leverett equations (2.1)-(2.3), then s does not reach the value 1. • For all θ in ∈]0, 1[, the "full saturation condition" is asymptotically satisfied, i.e.

lim ε→0 g(θ in ) 1 -f (θ in ) = θ in .

Proof

• The proof is straightforward from the fact that

x > g(x) 1 -f (x)
, ∀x ∈ [0, 1[ so that the "full saturation condition" cannot be satisfied.

• Recalling the properties of f and g (see also

lim ε→0 Q in f (θ in ) + v 0 h(0) g(θ in ) = Q in ,
which means that the "full saturation condition" is asymptotically satisfied.

Remark 5.3 From Propositions 5.1 and 5.2, it is impossible to get a (possible) stationary saturation with value 1 at any location of the domain, because the "full saturation condition" never holds (for non-pathological values of θ in ). Nevertheless, since it is asymptotically attained, we infer that considering small values of ε increases our hope to compare the Buckley-Leverett / Reynolds solution to the Elrod-Adams solution. Unfortunately, the behaviour of both fluxes f and g is pathological as ε tends to 0. Indeed, the graph of functions x → f (x) and x → g(x) clearly shows the presence of a boundary layer at x = 1 when ε tends to 0. Moreover, we have already mentioned that the numerical simulation of the Buckley-Leverett problem for small values of ε becomes difficult. Now let us enter into the details of the numerical comparison. Numerical data related to the simulation of the Buckley-Leverett equation are the following ones: ⊲ The mesh grid has 1400 elements.

⊲ The CFL condition is given by ∆t ∆x M = 0.9.

⊲ The numerical stationary solution is attained for a relative error in the discrete L 2 norm with value 10 -24 .

We study the behaviour of the bifluid in the two mentioned cases (see the Appendix at the end of the paper). 5.2.1 Case (i): the liquid phase is adhering to the lower (moving) surface Fig. 8 (resp. 9) gives the pressure (resp. saturation) distribution obtained with the Buckley-Leverett / Reynolds model along with the solution of the Elrod-Adams model. We can see that when ε decreases, the Buckley-Leverett saturation tends to the Elrod-Adams saturation; nevertheless, it never reaches the value 1, as previously mentionned. For ε = 10 -3 , which corresponds to physical situations in lubrication theory, the Buckley-Leverett saturation is very near to the Elrod-Adams one, up to the full saturation area. When considering the pressure, we can see that the peak of the Buckley-Leverett pressure never reaches the one of the Elrod-Adams pressure; this is due to the sensitivity to the coefficients in the generalized Reynolds equation. Indeed, the fact that the Buckley-Leverett saturation never reaches the value 1 prevents the generalized Reynolds pressure from approaching the Elrod-Adams pressure. Nevertheless, it is observed that non-positive pressures in the Buckley-Leverett / Reynolds model tend to vanish as ε tends to 0. Another interesting point is to consider that, using the Elrod-Adams saturation in the generalized Reynolds equation allows to exactly obtain the Elrod-Adams pressure, although the main differnece between the Elrod-Adams saturation and the Buckley-Leverett saturation lies in the (nearly) satured region (in the first case, the value 1 is reached although, in the second case, it never reaches this value); this observation shows how the "full saturation condition" is important. However, it is reasonable to say that computations for ε = 10 -3 allow to identify cavitated and satured areas.

(x, z) ∈ R 2 , 0 < x < 1, 0 < z < h(x)
with both models: the Buckley-Leverett / Reynolds model and the Elrod-Adams model.

• Computation of the velocity for the Buckley-Leverett / Reynolds model: Let us denote (U, V ) the velocity in the bifluid. The following relationship has been derived in [START_REF] Paoli | Asymptotic behavior of a two fluid flow in a thin domain: from Stokes equations to Buckley-Leverett equation and Reynolds law[END_REF]:

η(x, z) ∂ 2 U ∂z 2 (x, z) = dp dx (x), (5.2)
U (x, 0) = v 0 , (5.3)
U (x, h(x)) = 0, (5.4) with η(x, z) = µ l , if 0 < z < s(x)h(x), µ g , if s(x)h(x) < z < h(x).
Moreover, z → U (x, z) is continuous at the interface z = s(x)h(x). Thus, integrating twice Equation (5.2) and taking account of boundary conditions (5.3) and (5.4) allow to get the velocity profile, up to the velocity at the free boundary, denoted U ⋆ .

Nevertheless, U ⋆ is obtained by taking account of

lim z→(s(x)h(x)) -µ l ∂U ∂z (x, z) = lim z→(s(x)h(x)) + µ g ∂U ∂z (x, z) .
• Computation of the velocity for the Elrod-Adams model: The velocity profile for the Elrod-Adams model is computed under the additional assumption that the liquid phase of the lubricant is adhering to the moving surface (we recall that, actually, the Elrod-Adams model does not provide any information on the position of the liquid phase). Thus, the velocity at the free boundary (which separates the lubricant and the "empty" phase) is equal to 0. The velocity in the "empty" phase is reduced to 0. In details, we obtain:

µ l ∂ 2 U ∂z 2 (x, z) = dp dx (x), for z ∈]0, s(x)h(x)[
with the boundary conditions U (x, 0) = v 0 and U (x, s(x)h(x)) = 0. In the empty phase, we have U (x, z) = 0, for z ∈]s(x)h(x), h(x)[.

5.2.2

Case (ii): the liquid phase is adhering to the upper (fixed) surface Fig. 11 (resp. 12) gives the pressure (resp. saturation) distribution obtained with the Buckley-Leverett / Reynolds model for various values of ε. We can see that when ε decreases, the Buckley-Leverett saturation tends to the uniform distribution s ≡ θ in .

The pressure tends to 0 uniformly. Therefore, the Buckley-Leverett / Reynolds model does not approach the Elrod-Adams model in that configuration: assuming that the lubricant is adhering to the upper (fixed) surface is not relevant. In fact, we can observe that, in this configuration, the shear effects tend to vanish, as it is pointed out further. Fig. 13 describes the horizontal velocity distribution in the whole domain with the two models (Buckley-Leverett / Reynolds and Elrod-Adams).

• Computation of the velocity for the Buckley-Leverett / Reynolds model: Equations (5.2)-(5.4) still hold, up to this adapted definition (which takes into account the position of each phase):

η(x, z) = µ g , if 0 < z < (1 -s(x))h(x), µ l , if (1 -s(x))h(x) < z < h(x).
Now, as in Case (i), this is completed with the following equality at the interface:

lim z→((1-s(x))h(x)) - µ g ∂U ∂z (x, z) = lim z→((1-s(x))h(x)) + µ l ∂U ∂z (x, z) .
• Computation of the velocity for the Elrod-Adams model: The velocity profile for the Elrod-Adams model is computed under the additional assumption that the liquid phase of the lubricant is adhering to the upper surface (which is not supposed by the model).

We obtain:

µ l ∂ 2 U ∂z 2 (x, z) = dp dx (x), for z ∈](1 -s(x)) h(x), h(x)[,
with the boundary conditions U (x, h(x)) = v 0 and U (x, (1 -s(x)) h(x)) = 0. In the empty phase, we have U (x, z) = 0, for z ∈]0, (1 -s(x)) h(x)[.

Let us recall that the lubricant in liquid phase is considered as the reference fluid. Thus, velocity profiles show that the lubricant does not support the shear effects at all. This can be viewed also through the expression of the Buckley-Leverett flux, Q in f (s)+v 0 h(x)g(s), in which it is observed that when ε tends to 0, the non-classical contribution to the Buckley-Leverett flux, i.e. the shear contribution v 0 h(x)g(s), tends to vanish.

Conclusions

In this paper, a bifluid model (which has been derived under some geometrical assumption on the free boundary) has been studied in order to justify the Elrod-Adams model which is commonly used in tribology for the modelling of cavitation. In particular, we emphasize that the rupture of the continuous (liquid) film is recovered with the two approaches, which can be analyzed in terms of stationary entropy shock. Additionally, some differences still appear when comparing both solutions (pressures in particular), but this is due to little differences between the saturations linked to the coefficients sensitivity in the generalized Reynolds equation. However, the major obstacle which prevents us from formally identifying both models lies in the fact that the bifluid saturation never reaches the value 1 (although it is very close to this value). It can be thought that the Elrod-Adams model is attained as a singular limit (as the ratio of the viscosities ε tends to 0) of the bifluid model. This is of course an open question that we will focus on in forthcoming works.

Appendix. Expression of the flux functions

For the sake of simplicity, let us introduce the following function:

α (i) ε (s) = 1 -(1 -ε)s i , i = 1, 2, 3.
Case (i): the lubricant is adhering to the lower (moving) surface

This corresponds to Fig. 1. The functions A, B are given by: A(s) := A

(1)

ε (s) = 4α (1)
ε (s)α

ε (s) -3(α

ε (s)) 2 ε α

Figure 1 .

 1 Figure 1. Case (i). The liquid phase is adhering to the lower surface

Figure 2 .

 2 Figure 2. Case (ii). The liquid phase is adhering to the upper surface

Theorem 3 . 1 (

 31 Maximum principle) In Definition 1 and Equations (3.1)-(3.3), if we replace the functions f , g, Φ ±

  the flux is autonomous (for instance H ≡ 0) and if u 0 , u D are functions with values in an interval [a, b], then, u is a function with values in [a, b]; -For the auxiliary problem, if u 0 , u D are functions with values in an interval [a, b] ⊂ [0, 1], then, u is a function with values in [0, 1]; thus, only the set [0, 1] is preserved.

  for every ϕ ∈ H 1 0 (]0, T [×]0, 1[), so that Equation (3.1) is gained in a weak sense.

Theorem 3 . 6

 36 Under Assumption 5, problem (3.1)-(3.3) admits a unique weak entropy solution.

Remark 3 . 7 Y 8

 378 is an increasing function which defines an isomorphism from [0, L] to [0, 1]. In the same way, T is an increasing function which defines an isomorphism from [0, T ] to [0, T ], with Problem (2.1)-(2.3) can be reduced to an auxiliary problem (3.1)-(3.3).

Fig. 14 g

 14 and 15 in the Appendix), lim ε→0 = Id, uniformly on every compact set K ⊂ [0, 1[. Thus, we have (for θ in = 1):

Fig. 10

 10 Fig.10 describes the horizontal velocity distribution in the domain

,

  and the flux functions f and g are defined by: f (s) := f

  α

Figure 7 .Figure 8 .

 78 Figure 7. Saturation profile with partially active boundary conditions at different time steps: t = (i -1)T f /9, for i = 1 to 9, from left to right, top to bottom

Figure 9 .Figure 10 . 3 Figure 11 . 3 Figure 12 .Figure 13 . 4 Figure 14 .

 91031131213414 Figure 9. Saturation distribution for the Buckley-Leverett and Elrod-Adams models (Case (i))

Figure 15 . 4 Figure 16 . 4 Figure 17 . 4 Figure 18 . 4 Figure 19 .

 15416417418419 Figure 15. Shear contribution to the Buckley-Leverett flux function g (Case (i))

Acknowledgments

The authors are very grateful to Julien Vovelle for suggested ideas.

Case (ii): the lubricant is adhering to the upper (fixed) surface This corresponds to Fig. 2. We emphasize that the flux functions or Reynolds coefficients can be written, in a simple way, as a perturbation of the ones described in the first case. Indeed, the functions A, B are given by:

, and the flux functions f and g are defined by:

From the previous formulas, it can be deduced that only the shear terms are modified by the assumption on the geometrical assumption: indeed, in the Buckley-Leverett (resp. Reynolds) equation, only the flux function g (resp. right-hand side B) is modified.