A G Shvarts
email: andrei.shvarts@mines-paristech.fr

V A Yastrebov

Trapped fluid in contact interface

Keywords: trapped fluid, contact, surface roughness, local and global coefficient of friction, linear elastic fracture mechanics

Introduction

The study of mechanical contact and friction is a subject of high importance in many fields, from biological and engineering applications to geological sciences. Since natural and industrial surfaces always possess roughness under certain magnification, the contact between solid bodies occurs on separate patches corresponding to asperities of contacting surfaces, [START_REF] Archard | Contact and rubbing of flat surfaces[END_REF][START_REF] Archard | Elastic deformation and the laws of friction[END_REF][START_REF] Bowden | The friction and lubrication of solids[END_REF][START_REF] Greenwood | Contact of nominally flat surfaces[END_REF]. The evolution of the ratio of real contact area to apparent one under increasing external load determines essential contact properties such as friction, wear, adhesion, and is responsible for heat transport through contact interfaces. At the same time, the distribution of the free volume between contacting surfaces governs the fluid transport along the interface and is responsible for leakage/percolation phenomena, see for example [START_REF] Dapp | Self-affine elastic contacts: percolation and leakage[END_REF][START_REF] Paggi | Evolution of the free volume between rough surfaces in contact[END_REF]. Figure 1: Morphology of the contact interface between an elastic half-space with a rough surface and a rigid flat under increasing external load, numerical simulation results [START_REF] Yastrebov | From infinitesimal to full contact between rough surfaces: Evolution of the contact area[END_REF]: black is the real contact area, white is the "free" out-of-contact area and red is the "trapped" out-of-contact area, bounded inside non-simply connected contact patches.

Lubrication, i.e. separation of contacting surfaces by a fluid lubricant, is an efficient mechanism for friction and wear reducing. However, if the applied external load, pushing the contacting bodies together, is high enough or if the sliding velocities are small, the hydrodynamic pressure developing in the fluid is not sufficient to separate the solids, and asperities of both surfaces can get in direct contact despite the presence of the lubricant, which inevitably increases friction. This scenario corresponds to the so-called mixed regime, at which the load-bearing capacity is split between the fluid and the contact areas. For even higher pressures and lower velocities, the whole load is carried by the mechanical contacts, this regime is termed as boundary lubrication, see [START_REF] Hamrock | Fundamentals of fluid film lubrication[END_REF][START_REF] Azushima | Tribology in Sheet Rolling Technology[END_REF] for details. On the other hand, under increasing external load the lubricating fluid may be trapped in valleys (pools) delimited completely by the contact zone. Fig. 1 shows an example of the morphology of the contact interface between two elastic half-spaces with rough surfaces under external load [START_REF] Pei | Finite element modeling of elasto-plastic contact between rough surfaces[END_REF][START_REF] Carbone | Asperity contact theories: Do they predict linearity between contact area and load?[END_REF][START_REF] Putignano | The influence of the statistical properties of self-affine surfaces in elastic contacts: A numerical investigation[END_REF][START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF]. Note that the fraction of the "trapped" out-of-contact area (highlighted by red color), surrounded by contact patches, is significant.

The entrapment of the fluid in the interface can have a strong effect on the contact properties, especially if the fluid is highly incompressible [START_REF] Persson | Elastic contact mechanics: percolation of the contact area and fluid squeeze-out[END_REF][START_REF] Matsuda | Real contact area and friction property of rubber with two-dimensional regular wavy surface[END_REF]. First, the trapped fluid resists the compression, and thus opposes the growth of the real contact area. Second, the applied external load is shared between contacting asperities of the bodies and the pressurized fluid, so that the trapped fluid provides an additional load-carrying capacity (even in motionless contacts), reducing the normal pressure in the contact spots between the solid bodies. According to Coulomb's law of friction, the maximal tangential traction at the contact spots is proportional to the normal pressure, therefore the maximal macroscopic frictional force (of the whole contact interface) is proportional to the integral value of the normal pressure over the real contact area. Consequently, by taking into account the presence of the pressurized trapped fluid, a reduction of the global (apparent) coefficient of friction should take place.

The effect of lubricant entrapment on reduction of friction was first recognized in the study of cold metal forming processes [START_REF] Kudo | A note on the role of microscopically trapped lubricant at the tool-work interface[END_REF][START_REF] Nellemann | Real area of contact and friction stress -The role of trapped lubricant[END_REF], where the authors performed experiments on the sheet metal drawing test and identified three states, corresponding to different levels of the external pressure [START_REF] Azushima | Direct Observation of Contact Behaviour to Interpret the Pressure Dependence of the Coefficient of Friction in Sheet Metal Forming[END_REF]. Low values of external load are supported completely by the mechanical contact between asperities, and both global and local coefficients of friction are equal. At medium range of pressures, the global coefficient of friction decreases with increasing load due to closing of lubricant pools and generation of hydrostatic pressure in the fluid, which supports a part of the external load. At even higher load, fluid escapes from the pools and permeates into the contact zones, so that both the real contact area and the coefficient of friction decrease with increasing load. This effect is however biased by the fact that the real contact area does not evolve linearly under high pressures [START_REF] Archard | Elastic deformation and the laws of friction[END_REF], but rather as a concave function of pressure [START_REF] Persson | Theory of rubber friction and contact mechanics[END_REF][START_REF] Yastrebov | The role of the roughness spectral breadth in elastic contact of rough surfaces[END_REF], thus also resulting in formal decrease of the friction coefficient in contact spots. Nevertheless, experimental results together with finite-element simulations of the problem of entrapment and permeation of the fluid into the contact interface during upsetting of a cylinder were presented and aforementioned states were also identified [START_REF] Azushima | FEM Analysis of Hydrostatic Pressure Generated Within Lubricant Entrapped Into Pocket on Workpiece Surface in Upsetting Process[END_REF][START_REF] Azushima | Permeation of Lubricant Trapped Within Pocket Into Real Contact Area on the End Surface of Cylinder[END_REF]. An extensive experimental study of lubricant entrapping and escape in cold rolling processes was presented in [START_REF] Bech | Entrapment and escape of liquid lubricant in metal forming[END_REF].

In biological sciences the effect of trapped lubricant in human joints was investigated in the view of reduction of friction between rough cartilage surfaces [START_REF] Soltz | Hydrostatic Pressurization and Depletion of Trapped Lubricant Pool During Creep Contact of a Rippled Indenter Against a Biphasic Articular Cartilage Layer[END_REF][START_REF] Chan | The role of lubricant entrapment at biological interfaces: Reduction of friction and adhesion in articular cartilage[END_REF]. The concept of trapped fluid rises in the study of fatigue cracks in the rolling contact, which considers the process of crack growth due to pressurized fluid lubricant, forced inside of the crack by the external load and trapped there [START_REF] Bower | The Influence of Crack Face Friction and Trapped Fluid on Surface Initiated Rolling Contact Fatigue Cracks[END_REF]. The trapped fluid problem is also relevant to the geophysical studies: a landslide or an earthquake can be caused by an elevation of the pressure of the fluid in the pores inside the rock, see for example (Viesca and Rice, 2012;Garagash and Germanovich, 2012). The effect of the trapped fluid is also of interest for the study of basal sliding of glaciers, [START_REF] Cuffey | The physics of glaciers[END_REF]: the melt water, which is responsible for the lubrication, flows in a linked system of cavities in the interface between the glacier and the bedrock, and may be trapped there. Finally, the trapped fluid problem is also of importance for poromechanics [START_REF] Yu | A fractal permeability model for bi-dispersed porous media[END_REF][START_REF] Dormieux | Micromechanical approach to the behavior of poroelastic materials[END_REF][START_REF] Budiansky | Elastic moduli of a cracked solid[END_REF][START_REF] Coussy | Poromechanics[END_REF]. [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF] extended the Westergaard's celebrated analytical solution for the problem of contact between a regular wavy surface and a rigid half-plane (Westergaard, 1939) by taking into account the presence of a compressible fluid, trapped in the valleys between contacting asperities. Kuznetsov's solution demonstrates how the external pressure is divided between the fluid and the solid contact, which results in the decrease of the global coefficient of friction under increasing external load. However, due to the assumptions (i) that the wavy surface behaves as a flat one and that Flamant's solution holds for every surface point, and (ii) that the horizontal component of the fluid pressure is negligible, it cannot describe the escape of the lubricant and depletion of the real contact area. This question will be discussed later in detail. Recently an analytical solution was proposed for the problem of sliding of a rigid periodical punch along a viscoelastic Winkler's foundation with the incompressible fluid present in the gap [START_REF] Goryacheva | Modelling of a punch with a regular base relief sliding along a viscoelastic foundation with a liquid lubricant[END_REF].

Despite a significant attention to the problem of the trapped fluid in the contact interface, a few questions remain open, such as: the mechanism of the trap opening, the evolution of the real contact area and of the global coefficient of friction during this process, and also the distribution of the frictional shear tractions in the contact interface under external normal loading in the presence of the pressurized fluid in the interface. Note that these questions cannot be addressed in the framework of the boundary element method (BEM), since it assumes infinitesimal slopes of the surface roughness, which is, as we will show, a too restrictive assumption for the considered problem. We address these questions in the current study in the framework of the finite element method (FEM).

The paper is organized as follows. In Section 2 we present the statement of the problem of the mechanical contact coupled with pressurized compressible and incompressible fluids trapped in the contact interface. In Section 3 we outline existing analytical solutions of this problem, and in Section 4 we discuss methods for its numerical solution. Section 5 is devoted to results, including comparison of Kuznetsov's analytical solution with our numerical simulations, the evolution of the real contact area and the global coefficient of friction, as well as the simulation and analysis of the frictional behaviour of the system under normal and tangential external loading. In Section 6 we present the conclusions. Under the assumption of small slopes, in Appendix A we derive an analytical solution for vertical displacement of a wavy surface under the action of uniformly distributed pressure. Appendix B provides details of the numerical formulation of the coupled trapped fluid/mechanical contact problem.

Problem statement

We consider a mechanical contact problem between a deformable half-plane with a periodic wavy surface and a rigid flat under the action of a far-field external pressure. This case was historically the starting point for the study of contact of rough surfaces (Westergaard, 1939;[START_REF] Johnson | The contact of elastic regular wavy surfaces[END_REF]. In addition, we take into account the influence of compressible or incompressible fluid trapped in the free volume between the two bodies, see Fig. 2. We assume the plane strain problem and a linear elastic or elastic-perfectly J2-plastic (the latter one is presented in Section 5.4) isotropic constitutive laws for the solid. Note that this problem is similar to the one already solved by [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF], with the difference that we assume small but finite profile's slope, which, as will be shown in the paper, is of great importance for an accurate treatment of this problem. The initial gap between the wavy profile and the rigid plane, as well as the volume of this gap for one wavelength of the profile, are given, respectively, by:

g 0 (X) = ∆ 1 -cos 2πX λ , V g0 = l λ 0 ∆ 1 -cos 2πX λ dX = l λ ∆, (1)
where ∆ and λ are the amplitude and wavelength of the wavy surface profile, respectively, X is the horizontal coordinate in the initial (reference) configuration and l is the length in the direction of the third coordinate z, which under the assumption of the plane strain state of deformation will be assumed hereinafter equal to one length unit.

Analytical solutions

Westergaard's solution

The problem of contact between an elastic half-space with a regular wavy surface y = ∆ cos(2πx/λ) and a rigid flat without fluid in the interface was solved by Westergaard (1939), see also [START_REF] Johnson | The contact of elastic regular wavy surfaces[END_REF], for the case of infinitesimal ratio ∆/λ 1, i.e. infinitesimal slope of the roughness profile. According to this solution the pressure distribution inside contact patches (-a + λn ≤ x ≤ a + λn, n ∈ Z) is given by:

p W (x, a) = 2πE 1 -ν 2 ∆ λ cos πx λ sin 2 πa λ -sin 2 πx λ , ( 2
)
where a is the half-length of contact patch within one wavelength of the profile λ, and elsewhere p W = 0. E and ν are Young's modulus and Poisson's ratio, respectively. The mean pressure over the whole contact interface is given by pW

(a) = 1 λ λ 0 p W (x, a) dx = p * sin 2 πa λ , (3)
where p * = πE * ∆/λ is the pressure necessary to bring the entire interface in contact. In the static equilibrium pW is equal to the value of the external pressure that we will denote by p 0 . The complete contact is ensured, if p 0 ≥ p * . By introducing the notations A = 2a and A 0 = λ for the real and apparent contact areas, respectively, the ratio of the real contact area to the apparent one, based on the Westergaard's solution, is given by:

A A 0 = 2a λ = 2 π arcsin p 0 p * , 0 ≤ p 0 ≤ p * . (4)
3.2. Kuznetsov's solution [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF] extended the Westergaard's solution (2) by taking into account compressible fluid trapped in the valleys between contacting peaks of the wavy profile. Similarly, under the assumption of infinitesimal slope of the profile1 , the stress state in the contact interface in the presence of the additional fluid pressure, applied beyond the contact patches, was considered as the superposition of the stress state corresponding to the same contact area, but without influence of the fluid (i.e. the Westergaard's solution (2)), and a uniform field of the fluid pressure p f , applied everywhere and assumed not to distort the surface profile:

p K (x, a) = p f (a) + p W (x, a), if -a + λn ≤ x ≤ a + λn, n ∈ Z p f (a),
elsewhere.

(

Integration of p K (x, a) over one period of the waviness gives the following relation between the external pressure p 0 and the contact area: p 0 (a) = p f (a) + pW (a), where pW (a) was defined in (3).

The fluid pressure p f can be related to the current contact half-width a using a model of the compressible fluid with a bulk modulus K, which is defined as the ratio of infinitesimal pressure increase to the relative decrease of the volume:

K = -V f dp f dV f . (6)
In the linear compressibility model bulk modulus is a constant coefficient of proportionality between the relative change of volume of fluid and the fluid pressure [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF]:

p f = K 1 - V f V f 0 , (7)
where V f 0 is the volume of the fluid in unpressurized state and a smaller volume V f corresponds to the fluid pressure p f . However, the linear model of compressible fluid (7) does not provide satisfactory results for most of the fluids used in real-life lubrication problems, since a significant dependence of the compressibility modulus K of fluid on the pressure p f takes place [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF]. The simplest model, and yet quite precise for most of lubricating fluids, which takes into account this dependence is the compressibility linearly evolving with pressure [START_REF] Nellemann | Real area of contact and friction stress -The role of trapped lubricant[END_REF][START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF]:

K = K 0 + K 1 p f , (8)
where K 0 , K 1 > 0 are model parameters. The linear dependence (8), substituted into (6), upon integration results in the following non-linear relation between the fluid pressure and its volume:

p f = K 0 K 1 V f V f 0 -K1 -1 . (9)
Finally, it can be noted that the volume of the pressurized fluid V f is equal to the volume of the gap between the contacting surfaces V g , which can be found from the displacement field of the Westergaard's solution [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF] and related to the current contact half-width a:

V g (a) = V g0 1 -sin 2 πa λ 1 -ln sin 2 πa λ , (10)
where V g0 = l ∆ is the initial gap, i.e., corresponding to a = 0. We generalize original results [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF] and allow a partial filling of the initial gap by the fluid, so that V f 0 = θ V g0 , 0 < θ ≤ 1. Therefore, if the current gap volume is bigger than the initial fluid volume, V g > V f 0 , i.e. V g /V g0 > θ, then the fluid is not yet pressurized, and Westergaards solution is valid:

p 0 (a) = πE * ∆ λ sin 2 πa λ . If V g < V f 0 ,
or, equivalently, V g /V g0 < θ, the equation connecting the contact area and the external load has the following form in the case of linear compressible fluid:

p 0 (a) = πE * ∆ λ sin 2 πa λ + K θ θ -1 + sin 2 πa λ 1 -ln sin 2 πa λ , if V g /V g0 < θ, (11)
and in the case of non-linearly compressible fluid:

p 0 (a) = πE * ∆ λ sin 2 πa λ + K 0 K 1 θ K1 1 -sin 2 πa λ 1 -ln sin 2 πa λ -K1 -1 , if V g /V g0 < θ. ( 12
)
It is important to note also that Kuznetsov's solution even in the case of an arbitrary large modulus of compressibility of the fluid shows the growth of the contact patches under the increasing load. Furthermore, in the limit of incompressible fluid K → ∞ it gives a constant value of the real contact area, which can be found from the equation V g (a) = V f 0 . Consequently, Kuznetsov's solution, based on the assumption of infinitesimal slope of the profile, cannot predict depletion of the real contact area and escape of the fluid from the trap, which we demonstrate in following sections dropping out the assumption of infinitesimal slopes.

Numerical methods

Mechanical contact

In case of the unilateral contact between a deformable body and a rigid flat with an outer normal ν, the motion of the body is constrained, which can be formalized upon introduction of the normal gap function g -a signed distance from the points on the surface of the deformable body to the rigid plane: g > 0, when the point is separated from the plane, g < 0, when the point penetrates the plane (which is forbidden), g = 0, when the point is on the plane.

We will denote by Γ the potential contact zone (the whole surface), by Γ c ⊂ Γ the active contact zone, where the normal surface traction σ n must be negative in non-adhesive contact, and by Γ \ Γ c the inactive zone, which is out of contact. The constraints governing the frictionless unilateral contact problem are known as the Hertz-Signorini-Moreau conditions [START_REF] Wriggers | Computational Contact Mechanics[END_REF]:

g ≥ 0, σ n ≤ 0, σ n g = 0 at Γ ⇔ g = 0, σ n < 0, at Γ c g > 0, σ n = 0, at Γ \ Γ c . (13)
Therefore the considered problem is the constrained minimization problem for the potential energy of the mechanical system Π(u), where u is the displacement field. This problem can be solved using the Lagrange multipliers method [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF], with the Lagrangian functional defined as:

L(u, λ c ) = Π(u) + Γc λ c g(u) dΓ c , (14)
where λ c ≤ 0 is the Lagrange multiplier function, the values of which are equivalent to the normal traction in the contact zone.

In case of frictional contact, along with Hertz-Signorini-Moreau conditions (13), additional frictional constraints must be included in the problem, such as Coulomb's law of friction, which defines the following possible active contact states: • Slip:

• Stick: ġt = 0, |σ t | < µ |σ n |; g X 0 Ω 0 Γ ν ex ey n (a) p 0 Γc σ n g x Ω ν σ n ex ey n (b)
σ t = µ |σ n | ġt /| ġt |;
where ġt is the sliding velocity in the tangential plane between the corresponding points of the two surfaces, σ t is the tangential contact traction and µ is the coefficient of friction (CoF).

In order to include frictional constraints in the formulated above constrained minimization problem, special methods must be used, such as the penalty method (combined with the return mapping algorithm) or augmented Lagrangian method, for details see [START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF].

Trapped fluid constraints

Geometrical constraint for incompressible fluid

The area of the gap between the contacting surfaces V g in the presence of trapped incompressible fluid of volume V f must satisfy the following geometrical constraint:

V g ≥ V f = const, V g (X + u) = Γ f g(X + u) d Γ f , (15)
where Γ f = Γ \ Γ c and Γ f is the projection of Γ f on the rigid plane. The trapped fluid may fill completely or partially the gap between the contacting surfaces, therefore it can be present in two different states: "inactive", when V f < V g and the fluid if not pressurized (p f = 0), and "active", when V f = V g , and pressure in the fluid p f > 0, see Fig. 4a,4b. We may formulate this two states in a way similar to Hertz-Signorini-Moreau conditions: 16)

V g ≥ V f , p f ≥ 0, p f (V g -V f ) = 0 ⇔ V g = V f , p f > 0, (active state) V g > V f , p f = 0, (inactive state). (

Simulation of incompressible fluid using a Lagrange multiplier

In the case of the inactive state of the trapped fluid we have only the mechanical contact problem between the elastic body and the rigid plane, while if the fluid is in the active state, we must consider additionally the gap volume constraint (15). The Lagrange multiplier method may be used again in order to fulfill this constraint, and the combined functional for the coupled problem can be defined as:

L(u, λ c , λ f ) = Π(u) + Γc λ c g(u) dΓ c -λ f (V g (u) -V f ), ( 17
)
where λ f ≥ 0 is the Lagrange multiplier for the trapped fluid problem, which is equivalent to the fluid pressure p f . The solution of the coupled problem is a stationary point of the Lagrangian (17), which requires the calculation of its variation:

δL(u, λ c , λ f ) = ∂Π(u) ∂u δu + Γc δλ c g(u) + λ c ∂g(u) ∂u δu dΓ c -δλ f (V g (u) -V f ) + λ f ∂V g (u) ∂u δu = 0. ( 18
) V f Ω Γ Γ f V g >V f ex ey (a) p 0 p f σ n Ω ex ey σ n V g =V f (b)

Simulation of the compressible fluid with the penalty method

The geometrical constrain (15) for the trapped fluid can also be treated with the penalty method. In accordance with the linear penalty method, instead of the term λ f (V g (u) -V f ), the following term should be added in ( 17) to take into account the trapped fluid constraint:

W f (u) = 2 (V f 0 -V g (u)) 2 , ( 19
)
if the fluid is in active state V g < V f 0 , and zero otherwise. In the above formula is the penalty parameter, and V f 0 is the initial volume of the fluid. Let us assume that the fluid is in active state. Calculating the variation of ( 19), we obtain the contribution of the trapped fluid to the balance of virtual works:

δW f (u) = -(V f 0 -V g (u)) ∂V g (u) ∂u δu, (20)
where the value of the term (V f 0 -V g (u)) equals to the fluid pressure p f . Under the penalty formulation the gap volume constraint ( 15) is never satisfied exactly, i.e. the current volume of the active fluid V g (u) is always smaller, than the initial fluid volume V f 0 . Therefore, the penalty method corresponds to the model of the compressible fluid, and a comparison between ( 7) and ( 20) shows that the linear penalty method represents the compressibility model with the constant bulk modulus K, if = K/V f 0 .

In order to simulate the behaviour of the compressible fluid with pressure-dependent bulk modulus ( 8)-( 9), the non-linear penalty method for the trapped fluid constraint (15) may be used. The contribution of the fluid to the balance of virtual works in this case takes the form:

δW f = - K 0 K 1 V g (u) V f 0 -K1 -1 ∂V g (u) ∂u δu. (21)

Results and discussion

We solved the coupled problem using the finite element method with implemented monolithic coupling scheme in finite element suite Z-set [START_REF] Besson | Large scale object-oriented finite element code design[END_REF][START_REF]Nonlinear material and structure analysis suite[END_REF]. Contrary to Kuznetsov's analytical results or BEM analyses, we did not assume infinitesimal slopes, i.e. the value ∆/λ is arbitrary. We used a finite element mesh with 1024 nodes in the contact interface per wavelength (19364 nodes in total in the structural mesh), see Fig. 5. Hereinafter, if not mentioned differently, we considered the roughness profile with ∆/λ = 0.01. In the following, we will also discuss how this ratio affects the results. The horizontal dimension of the finite element mesh equals to the wavelength λ and the ratio of the profile amplitude ∆ to the vertical mesh dimension H is ∆/H = 0.005. On the vertical boundaries of the mesh we apply symmetry boundary conditions (u x = 0), and the bottom edge of the deformable solid is displaced vertically towards the rigid flat within 200 load steps. A corotational updated Lagrangian framework was used in our simulations, which is needed to capture properly that the fluid pressure applied to the updated configuration is collinear to element normals. In simulation we measure the vertical reaction, the extension of the contact area, the pressure in the contact zone and the fluid pressure. Hereinafter, if not mentioned differently, we performed frictionless simulations, and estimated the value of the global coefficient of friction using the following approach. We consider the global coefficient of friction as the coefficient of proportionality between the maximal tangential force per wavelength F t and the normal one We neglect shear forces in the fluid and therefore the ratio between the global and local coefficients of friction can be calculated as:

F n = p 0 λ, i.e. |F t | ≤ µ glob |F n |.
µ glob µ loc = Γc |σ n | dΓ c |F n | = Γc |σ n | dΓ c p 0 λ. (22)
Finally, using notations for the real A and the apparent A 0 contact areas, which were introduced above, Eq. ( 22) can be rewritten as:

µ glob µ loc = 1 - p f p 0 1 - A A 0 . (23)

Incompressible fluid

In this section we study the model of an incompressible fluid trapped in the contact interface. Note that real-life lubricating fluids have significantly lower initial bulk moduli than metals. Nevertheless, this idealized model enables us to focus on the mechanism of the trap opening by the pressurized fluid, while compressible fluids will be considered in the following sections.

V f / V g0 0.7 0.8 0.9 0.7 0.8 0.9 0.7 0.8 0.9 We study the evolution of the real contact area in the presence of incompressible fluid in the interface under the increasing external pressure using the Lagrange multiplier method. We investigate how the magnitude of the slope of the profile (∆/λ) and the ratio between the trapped fluid volume and the initial gap volume V f /V g0 affect the solution of the coupled problem. The distribution of some stress and strain components in the bulk of the deformable solid during the process of trap opening is shown in Fig. 8.

The evolution of the contact area close to the moment of the activation of the fluid is presented in Fig. 6(a). The regime in which the fluid is not yet pressurized (V g > V f ) coincides with Westergaard's equation (2). According to this analytical solution the ratio of the current volume of the gap to the initial one V g /V g0 is a monotonically decreasing function of contact area and does not depend on the slope of the profile ∆/λ, see (10). Therefore, the contact area A act , reached when the fluid gets pressurized (V g = V f ) does not depend on the slope of the profile and is increasing with decreasing V f /V g0 . For a given ∆/λ the pressure necessary to activate the fluid p act is also increasing with decreasing V f /V g0 . At the same time, for a given V f /V g0 , the value of p act is proportional to the slope ∆/λ.

One can note in Fig. 6(a) that once the fluid is pressurized, the contact area is slowly decreasing, contrary to the Kuznetsov's solution, which predicts the contact area to remain constant. In Fig. 6(b) we show the evolution of the contact area in a much wider range of loads, than in Fig. 6(a), and observe a monotonic decrease of the contact area, ultimately it reaches zero value, which corresponds to the opening of the trap. Surprisingly, results of simulations with different (decreasing) profile slope ∆/λ do not tend to the Kuznetsov's solution (derived under assumption of infinitesimal ∆/λ and assuming that the wave profile is similar to a flat one), but converge to a different limit! At the same time we observe that the external pressure necessary to open the trap p open also converges to a certain limit with ∆/λ → 0.

In order to explain this intriguing result, first we note that since the solution of linearly elastic problem with and without contact is unique, see, for example, [START_REF] Truesdell | Mechanics of Solids: Volume II, Linear Theories of Elasticity and Thermoelasticity, Linear and Nonlinear Theories of Rods, Plates, and Shells[END_REF] Evolution of the ratio V f /V g0 (V f is the volume between the deformed surface and a horizontal plane y = y 0 , where y 0 is the current position of the crest, and V g0 is the initial volume of the gap) with the increasing hydrostatic pressure p f for several profiles with different slope ∆/λ.

an auxiliary problem of the uniform hydrostatic fluid pressure on the wavy profile, see Fig. 7(a). The Kuznetsov's solution is based on an assumption that a uniform distribution of the hydrostatic pressure does not distort the wavy surface [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF]. In our numerical simulations we showed that for small, but finite ∆/λ this assumption does not hold, the wavy surface distorts: the crest's displacement is bigger than the displacement of the trough, which is quite an evident result.

Due to the non-zero slope of the contact interface, the fluid pressure acts not only in the vertical direction but also in the horizontal one, thus leading to the additional in-plane compression of the material near the crest and, on the opposite, to the additional in-plane tensile contribution near the trough, see Fig. 7(a). Thus, there exists a linearly elastic solution for a uniformly distributed pressure p f , which results in such surface deformation, that the integral of the gap equals to the fluid volume V f , i.e.:

∃p f such that Γ (y 0 -(X y + u y ))dΓ = V f , (24)
where y 0 is the position of the crest after applying the uniform pressure p f . We derived an analytical formula for computation of V f , based on the assumption of small, but finite ∆/λ:

V f /V g0 = 1 - 2(1 -2ν)(1 + ν)p f E , (25)
see Appendix A for details. The relative change of volume induced by a uniformly applied pressure p f does not depend on the value ∆/λ, but only on elastic properties of the solid. In Fig. 7(b) the comparison of this formula with the numerical results for several profiles with different ∆/λ is shown. Numerical results are tending towards the analytical solution with decreasing ∆/λ. Therefore, we have shown that for any given V f /V g0 there exists uniform pressure p f , which results in a such distortion of the surface, that the volume between the surface and a rigid flat equals to V f . Moreover, in the limit of infinitesimal slopes, this critical pressure does not depend on the slope. Equation ( 25) can be readily used to compute the pressure needed to open the trap; it is valid for both incompressible and compressible fluids.

The obtained result explains why the curves of evolution of the real contact area with the increasing pressure for surfaces with different slopes tend to a certain limit with decreasing ∆/λ (which remains, -7e+03 -6.6e+03 -7e+03 -6.9e+03 -6.9e+03 -6.9e+03 -6.8e+03 -6.8e+03 -6.7e+03 -6.7e+03 -6.7e+03 -6.6e+03 -6.6e+03 -6.9e+03 -6.9e+03 -6.9e+03 -6.8e+03 -6.8e+03 -6.7e+03 -6.7e+03 -6.7e+03 -6.6e+03 -6.6e+03 sig22 map:67.0000 time:0.67 min:-7006.45 max:-6564.98 -6.3% ε yy -9.3% -8.6% ε yy -5.7% V f / V g0 =0.9

-1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 sig22
A/A 0 =9% A/A 0 =4.5% A/A 0 =0% -0.7% 0% -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 -1e+04 sig22
V f / V g0 =0.8 V f / V g0 =0.7 (a) (b) V f / V g0 =0.9 V f / V g0 =0.8 V f / V g0 =0.7 A act /A 0 (V f /V g0 =0.8) A act /A 0 A act /A 0 (V f /V g0 =0.7) V f / V g0 =0.9 V f / V g0 =0.8 V f / V g0 =0.7 (c) (d) Kuznetsov's approximation (V f /V g0 =0.7)

Frictionless simulation

Frictional simulation Kuznetsov's solution however, finite) see Fig. 6(b), while p open is different for different V f /V g0 , see Fig. 9(a). In the case when the fluid fills completely the initial gap volume (not shown in figures), i.e. V f = V g0 , the real contact area equals zero during the whole process of loading, i.e. the fluid is pressurized, but is never trapped. We present in Fig. 9(b) the evolution of the integral of normal tractions over the contact area, normalized by λ and E * , for different values of ∆/λ and V f /V g0 . The results show that just after the fluid becomes pressurized, the integral of contact pressure grows almost linearly, which is also predicted by the Kuznetsov's solution in the limit K → ∞:

1 E * λ Γc |σ n | dΓ c = p 0 E * A act A 0 + π 1 - A act A 0 ∆ λ sin 2 π 2 A act A 0 , (26)
where, contrary to numerical results, it was assumed that A act remains constant under the increasing external pressure p 0 . However, due to the fact that we consider finite slope of the profile in the numerical solution, the linear part in the dependence of contact reaction on external pressure is followed by a non-linear concave part, reaching maximum value and then decreasing to zero. Consequently, the global coefficient of friction also vanishes. The results on the estimation of the ratio between global and local coefficients of friction are presented in Fig. 9(c). Before the fluid gets pressurized, the global CoF equals to the local one. After that, the global CoF is monotonically decreasing with the increasing external pressure p 0 . This decrease is related to repartition of the external load between the contact and the fluid; the latter is assumed not to resist shear in the quasi-static limit. Note that for high values of p 0 , i.e. close to opening of the trap, the evolution of the global CoF is independent from the slope (∆/λ) and depends only on the ratio V f /V g0 . On the other hand, for low values of p 0 slightly higher than the activation pressure (see Fig. 9(d)) the Kuznetsov's solution under the assumption of infinite K shows the global CoF decreasing as 1/p 0 :

µ glob µ loc = A act A 0 + π 1 - A act A 0 ∆ λ E * p 0 sin 2 π 2 A act A 0 . (27)
Note, that the term containing 1/p 0 is proportional to the ratio ∆/λ. In addition to estimations of the global coefficient of friction ( 22)-( 23), based on the frictionless simulation of the coupled problem under normal loading, we performed the direct computation of µ glob = |F t |/|F n | in the frictional simulation of the coupled problem during sliding under normal and tangential loads. Note that in the latter simulation for both normal and frictional contact constraints we use the augmented Lagrangian method and the classic Lagrange multiplier method for the fluid constraint. The comparison of the results is presented in Fig. 9(d) for the case of V f /V g0 = 0.9 and different ratios of ∆/λ: Kuznetsov's solution ( 27) is presented with dashed curves, estimations based on frictionless simulation are shown as solid curves; results, which were obtained taking explicitly into account friction in the interface, are presented as crosses for a few particular values of external pressure p 0 . This comparison shows that the frictionless result, based on the assumption of separate consideration of tangential and normal contributions in the interface [START_REF] Johnson | Contact mechanics[END_REF], provides a trustworthy estimation of the global coefficient of friction.

Note that these considerations can be applied to multi-cracked materials such as rocks with fluid in contact interfaces. The irreversible deformation in rocks is related to the frictional sliding at crack interfaces, which starts after the mean shear traction σ t in the interface reaches the frictional limit determined by the coefficient of friction and the contact pressure µ glob σ n . Being homogenized over all randomly oriented crack orientations, these considerations give rise to Drucker-Prager-type constitutive behaviour with the initial yield surface given by f = σ vm + µ glob p -R 0 , where σ vm is the von Mises stress, p = -trace(σ)/3 is the hydrostatic pressure and R 0 is the initial yield stress for pure shear. Because of the presence of an incompressible fluid in the interface, the frictional limit does not increase linearly (or equivalently the global coefficient of friction does not remain constant), but reaches its maximum and decreases down to zero as shown in Fig. 9(b). This behaviour is very similar to advanced pressure-dependent plasticity models with a so-called cap, which corresponds to the decay of the von Mises yield stress with increasing pressure [START_REF] Resende | Formulation of Drucker-Prager cap model[END_REF]. But contrary to the pore-collapse mechanism [START_REF] Suarez-Rivera | Indentation by pore collapse in porous rocks[END_REF][START_REF] Perrin | Rudnicki and Rice's analysis of strain localization revisited[END_REF][START_REF] Issen | Conditions for compaction bands in porous rock[END_REF], here this decay results from the decrease of the global friction with the hydrostatic pressure in presence of the fluid, this result also holds for non-linearly compressible fluids.

Compressible fluid with constant bulk modulus

Here our analysis is extended to the case of compressible fluids. In Fig. 10(a) we present the comparison of the numerical simulation of a linearly compressible trapped fluid under the linear penalty formulation with the analytical solution (11). We plot the evolution of the ratio of the real contact area to the apparent one under increasing external pressure for the case when the fluid occupies 70% of the initial gap, i.e. V f 0 /V g0 = 0.7. Different curves correspond to different values of the modulus of compressibility of the fluid K f , normalized by the bulk modulus of the solid body K s = E/3(1 -2ν), and for each numerical result a corresponding analytical curve is presented for comparison. Before pressurization of the fluid, the presence of the latter does not affect the solution and all curves follow the Westergaard's solution (2). For the pressurized fluid, the results show a good agreement between numerical and analytical solutions for values K f /K s 1, and for K f ≈ 0 the solution coincides completely with the Westergaard's formula. However, with the increase of the K f , in the region corresponding to the active fluid, the difference between numerical and analytical solutions becomes more pronounced. For the ratio K f /K s close to unity, the numerical results shows an almost constant value of the real contact area under the increasing load. Note, that the same result will hold for an incompressible fluid trapped the interface between two incompressible solids.

K f /K s =0.001 K f /K s =0.05 K f /K s =0.1 K f /K s =1.0 K f /K s =100 (a) (b)
For even greater K f /K s , the numerical results show a decrease of the real contact area, which means that the pressurized fluid starts to open the contact. Due to inherent assumptions of infinitesimal slopes, these effects cannot be predicted by the analytical solution.

In Fig. 10(a) the results were presented for V f 0 /V g0 = 0.7, note that the smaller this ratio is, the bigger are the value of pressure necessary to bring the fluid in active state and the corresponding value of the contact area. However, after the fluid becomes pressurized, for sufficiently high values of external pressure, the evolution of the contact area is influenced only by the compressibility modulus of the fluid and the mean slope of the profile. The bigger is the compressibility modulus or the slope, the smaller is the contact area for the same external pressure.

To emphasize the difference between the analytical and numerical solutions for a nearly incompressible fluid, we plot the pressure distribution near a contact patch under the increasing load for both solutions, see Fig. 10(b). The representation of the stress state in the contact patches as a superposition of the stress state for the same contact area without the influence of the fluid a uniform fluid pressure, applied everywhere ( 5) still holds for the numerical solution, but unlike the analytic solution, in our results a significant reduction of the contact area for nearly incompressible fluid is observed. Note that in our numerical solution for sufficiently high external pressure the real contact area vanishes, which means that the fluid separates the contacting surfaces everywhere, and the external pressure is entirely supported by the fluid under the pressure equal to the external one p f = p 0 .

Compressible fluid with pressure-dependent bulk modulus

As was shown in Fig. 10(a) for the case of linearly compressible fluid (with constant bulk modulus), starting from the pressurization of the fluid, the real contact area evolves monotonically with the external pressure: if the fluid bulk modulus is considerably less than the one of the solid (K f K s ), then the real contact area monotonically increases, if K f K S , then the contact area decreases down to zero, corresponding to the opening of the trap. The latter case is interesting for the study of the process of the fluid permeation into the contact zone and reduction of the global coefficient of friction, however, as it was mentioned in the Sec. 5.1 for the incompressible fluid, the situation when the initial fluid bulk modulus is greater than that of the solid remains non-physical and serves as an idealized model. On the other hand, real fluids behave non-linearly and their bulk modulus increases with increasing pressure, and thus even if the fluid bulk modulus is smaller than that of the solid in the first stage of pressurization, it eventually becomes greater than the one of the solid under the increasing pressure.

We present results of the numerical simulation for coupled problem with non-linear fluids: evolution of the contact area and global coefficient of friction with increasing external pressure, see Figs. 11(a),(b), respectively. Physically relevant values for two solid materials are used: a typical steel (E = 200 GPa, ν = 0.28, K s ≈ 151.5 GPa) and aluminium (E = 70 GPa, ν = 0.33, K s ≈ 83.33 GPa), and three types of fluid (see Eq. ( 8)): water (K 0 = 2112.5 MPa, K 1 = 6.5), glycerine (K 0 = 4151.5 MPa, K 1 = 8.74) and a typical mineral oil (K 0 = 2000.0 MPa, K 1 = 9.25) [START_REF] Kuznetsov | Effect of fluid lubricant on the contact characteristics of rough elastic bodies in compression[END_REF][START_REF] Nellemann | Real area of contact and friction stress -The role of trapped lubricant[END_REF]. We limit this study to the contact problem with the fluid completely filling up the gap (but only up to the upper boundary) during the whole process of loading. Such formulation remains rather general since, due to the realistic fluid model, the contact zone will inevitably appear in the beginning of loading. At low external pressures numerical results coincide with the analytical solutions for non-linear fluids, see (12). However, in contrast to the analytical solution, which cannot account for depletion of the contact zone, the numerically obtained contact area reaches a certain maximal value and then decreases (the latter is observed for external pressures p 0 /E * > 0.2, not shown in figures since the contact area does not reach zero value even for extremely high values of the external pressure p 0 = E * ). Note that for each of considered materials, the obtained curves for water and oil coincide in the beginning of loading due to almost equal initial bulk moduli K 0 of these fluids, and deviate for higher external pressures due to difference in K 1 , while for glycerine K 0 is significantly bigger, leading to a smaller contact area in this case.

The global coefficient of friction (CoF) also shows a non-monotonic behaviour, see Figs. 11(b), first, it rapidly increases up to a certain maximal value. Within this stage, the numerical and analytical results are very close, while for higher pressures a strong deviation of analytical and numerical results is observed. In analytical solution, even though the global CoF may decrease after the first extremum-maximum (see results obtained for the steel), it eventually increases again after reaching the second extremum-minimum. More accurate numerical results predict a monotonic decrease of the global CoF after reaching the first maximum. Note that in the simplified case considered here, the effect of hydrostatic lubrication lowers significantly the maximal global CoF, which does not exceed ≈ 36 % of the local CoF for the steel, and does not exceed ≈ 24 % of the local CoF for the aluminium. Such a strongly non-linear behaviour of the global coefficient of friction is explained by a competition between non-linear fluid pressurization and non-linear contact area evolution (see Eq. ( 23)).

The numerical solution shows that the maximal value of the CoF and its slope after passing the extremum both depend on the ratio between the bulk moduli of the fluid K f = K 0 + K 1 p f and the solid K s . The bigger is the initial modulus K 0 , the higher is the maximal CoF (which explains almost equal peak values of the CoF for water and oil and much lower value for glycerine). At the same time, the bigger is the coefficient K 1 , the faster the CoF decreases.

Steel

Alum.

Steel

Alum. 12). Evolution of (c) the ratio of real contact area to the apparent one A/A 0 , and of (d) the global to local coefficients of friction under increasing external pressure in the case of elastic-perfectly plastic solid and incompressible fluid. Note that in the initial configuration the fluid does not occupy the entire gap: V f /V g0 = 0.9. Dashed curves are presented for comparison with the cases of purely elastic solids, discussed in Sec. 5.1. Vertical dash-dotted line indicates the hardness taken to be H = 3σ Y .

3σ Y /E * (steel) 3σ Y /E * (alum.) (d)
We performed additional simulations varying the slope of the roughness profile ∆/λ in the interval [0.005; 0.02]. The results showed that the evolution of the real contact area is almost independent of the ratio ∆/λ (similarly to the case of the incompressible fluid). On the other hand, variation of this ratio has a considerable effect on the peak value of the global CoF, which increases with increasing ∆/λ. However, for high values of external pressure, the CoF does not depend on the slope of the profile, as it was also observed for the incompressible case.

5.4. Elastic-perfectly plastic solid 0 1.7%

-7e+03 -6.6e+03 -7e+03 -6.9e+03 -6.9e+03 -6.9e+03 -6.8e+03 -6.8e+03 -6.7e+03 -6.7e+03 -6.7e+03 -6.6e+03 -6.6e+03 As was shown in the previous section, the pressurized fluid can open the contact only at unrealistically high pressures, at which the linearly elastic material model is irrelevant. Thus, here a more realistic case is presented taking into account a non-linearly compressible fluid and an elasto-plastic material behaviour. Moreover, such a study is relevant for the aforementioned experimental observations of entrapment and consequent permeation of the fluid into the contact interface in metal forming processes and during upsetting of an elasto-plastic cylinder [START_REF] Azushima | Direct Observation of Contact Behaviour to Interpret the Pressure Dependence of the Coefficient of Friction in Sheet Metal Forming[END_REF][START_REF] Azushima | Permeation of Lubricant Trapped Within Pocket Into Real Contact Area on the End Surface of Cylinder[END_REF]. Here we consider elasticperfectly plastic materials (von Mises stress criterion): steel, E = 200 GPa, ν = 0.28, yield stress σ Y = 250 MPa and aluminium, E = 70 GPa, ν = 0.33, σ Y = 240 MPa.

It is well known that in elasto-plastic mechanical contact, the contact pressure cannot exceed the material hardness, which can be reliably estimated as H ≈ 3σ Y [START_REF] Bowden | The friction and lubrication of solids[END_REF][START_REF] Johnson | Contact mechanics[END_REF][START_REF] Mesarovic | Spherical indentation of elastic-plastic solids[END_REF]. Thus it could be expected that after the pressure in the fluid reaches material hardness the contact abruptly opens. However, as demonstrated by our simulations, due to the high hydrostatic compressive state, the pressure in the contact can significantly overpass the material hardness.

First, we study incompressible fluid, and present in Fig. 11(c) the evolution of contact area in the case of V f /V g0 = 0.9. It shows significantly different behaviour compared to elastic material: after the fluid becomes activated, the contact area is non-monotonic function of external pressure, it has a small increase, and then an abrupt decrease, corresponding to the state when fluid pressure reaches the value of contact pressure, and, consequently, permeation becomes possible. Normal tractions in contact interface increase beyond 6σ Y due to hydrostatic pressurization of the solid. In Fig. 11(d) the resulting evolution of the global CoF is presented, which shows considerably lower values of the CoF for both considered materials, than the ones observed in the purely elastic case (for the same external pressure). Fields of the accumulated plastic strain in the solid at different loading steps are presented in Fig. 12, note that once the fluid gets pressurized, the plastic zone is not limited to the contact vicinity, but spreads over the entire interface and, consequently, the whole bulk of the solid. Notably, a secondary onset of plastic deformation appears in the trough of the wavy profile, it complements the classical plastic core appearing under the contact zone and spreading to the contact interface [START_REF] Johnson | Contact mechanics[END_REF][START_REF] Mesarovic | Spherical indentation of elastic-plastic solids[END_REF][START_REF] Kogut | Elastic-plastic contact analysis of a sphere and a rigid flat[END_REF]Alcalá and Esqué-de los Ojos, 2010).

Varying the slope of the profile as in Sections 5.1 and 5.3, we showed that in contrast to the case of elastic solids, where the evolution of the contact area during the process of trap opening does not depend on the slope of the profile ∆/λ, in case of elasto-plastic solids, for a given ratio V f /V g0 , once the fluid gets pressurized, the higher is the ratio ∆/λ, the bigger is the contact area.

The behaviour of the system incorporating the elasto-plastic material and non-linearly compressible fluid is shown in Figs. 13(a-d): the contact area after reaching its maximum abruptly decreases, resulting in a fast permeation of the fluid in the contact interface and eventual opening of the contact. Note that after a relatively fast saturation of the contact pressure at approximate material hardness H ≈ 3σ Y , a further increase in pressure without fluid permeation still remains possible up to huge pressure values p 0 σ y . In reality however, due to the micro-roughness permeation of the fluid in the contact interface may happen on earlier stages of the deformation.

In Fig. 13(c,d) the evolution of the global CoF is depicted, which shows a rather similar behaviour to the one observed in the case of the elastic solid, having multiple extrema in the beginning of loading. Note that the amplitude of the first maximum of CoF is increasing with increasing slope of the profile, which was also observed in the simulations with the purely elastic material.

Friction in the contact interface

Finally, we discuss the distribution of frictional tractions in the contact interface during the process of the trap opening. We consider a coupled problem for an incompressible fluid with Coulomb's friction in the contact interface, as in previous analysis the shear forces in the trapped fluid are neglected due to quasi-static analysis. The following geometrical parameters are used: ∆/λ = 0.01, V f /V g0 = 0.95. In order to obtain more reliable results, we refined the mesh to have 512 nodes within the maximal extension of the contact zone a/λ = 0.05, with 1024 surface elements in total.

We emphasize that the external tangential load is absent in this study, thus we consider here only normal loading, in which two stages can be distinguished. During the first stage the external pressure p 0 increases from zero value to p act , the value necessary to bring the fluid into active state, and the contact area reaches the maximum value. Results for the first stage are presented in Fig. 14(a), where, in order to visualize stick and slip zones, we plot normal tractions, multiplied by the coefficient of friction (CoF) µ = 0.2. Those results are very close to the classic self-similar (remaining the same for any load under a proper coordinate/pressure scaling [START_REF] Spence | Self similar solutions to adhesive contact problems with incremental loading[END_REF])) distribution of tractions, because the wavy profile in the region of interest is very close to a parabolic curve. During the second stage of loading (p 0 > p act ) the fluid is in the pressurized state and influences the interfacial traction distribution.

Since the slope of the roughness profile is small, the distribution of normal traction should resemble, at least for p 0 not much greater than p act , the analytical solution for a fluid bulk modulus tending to infinity (K → ∞), in which a uniform pressure offset is added everywhere to the field of the normal traction corresponding to the external pressure p act . In accordance to that, tangential traction remains almost unchanged over the majority of the contact interface. Since the contact pressure is increased by the fluid pressure offset, all points pass to the stick state, i.e. adhere to their positions. However, due to the finiteness of the slope being taken into account, the distribution of normal traction slightly differs from the analytical solution in the same way as was discussed in Sec. 5.2,see Fig. 10(b), i.e. a slight decrease of the contact area takes place.

For p 0 sufficiently greater than p act , see Fig. 14(b), the effects of finite slope become more pronounced, the contact area is gradually decreasing and a remarkable evolution of the tangential traction is observed. A singularity in the tangential traction emerges at the boundary of the contact zone, with the value at the tip of this singularity limited by the Coulomb's law. In order to explain and verify this intriguing result, we consider an analogy between the process of the trap opening with the interfacial friction and the mode-II crack propagation in the framework of linear elastic fracture mechanics (LEFM) theory (Tada et al., 1973).

Note that the analogy is not complete in physical sense: during the process of trap opening due to pressurization of the incompressible fluid, new surface is not created, since no atomic bonds must be broken in order to separate the surfaces. The physical reason for the singularity in tangential stress is the following: when points of the surface loose contact, their normal traction reduces not down to zero, but to the value of fluid pressure, thus the frictional limit near the contact edge remains elevated. Thus, the points of the interface before loosing contact have non-zero shear traction, and being liberated from this traction after loosing the contact, these points slide freely, in absence of frictional resistance, towards the centre of the contact zone.

The fluid activation corresponds to the maximal extension of the contact zone (we shall denote the maximal contact half-length as a * , and during the subsequent increase of the external pressure the width of the contact zone is monotonically decreasing. For sufficiently small slope of the roughness profile, the situation corresponding to contact half-length a < a * can be considered as a configuration of two bonded dissimilar solids with two aligned semi-infinite interfacial cracks in the interface, separated by 2a, see Fig. 14(c). Using the superposition principle, the observed stress state, corresponding to the half-length of the contact patch a, can be represented as a superposition of the initial shear traction σ * t (x), corresponding to the moment of activation of the fluid, and a stress induced by the same traction with the opposite sign, σ - t (x) = -σ * t (x) applied only on the surfaces of the cracks in the intervals x ∈ [-a * , -a] and [a, a * ]. Such traction induces a singular shear stresses in the region between two cracks x ∈ [-a, a], thus σ - t (x) can be written as:

σ - t (x) =        -σ * t (x), x ∈ [-a * , -a] ∪ [a, a * ] 1 √ 2π Im K(a, σ * t ) (x-a) i √ |x-a| -(x+a) i √ |x+a| , x ∈ [-a, a] 0, |x| > a * , (28)
where K is the complex stress intensity factor, see (Rice and Sih, 1965;[START_REF] Rice | Elastic fracture mechanics concepts for interfacial cracks[END_REF], and two terms in brackets in (28.2) correspond to two semi-infinite cracks being considered, so that σ - t (0) = σ * t (0) = 0, Im is the imaginary part. Therefore, the resulting distribution of shear tractions is given by the superposition σ t (x) = σ * t (x) + σ - t (x). The complex stress intensity factor K is calculated using the existing analytical formula for considered configuration and shear traction distribution (Rice and Sih, 1965;[START_REF] Rice | Elastic fracture mechanics concepts for interfacial cracks[END_REF]:

K(a, σ * t ) = [k 1 (a, σ * t ) + ik 2 (a, σ * t )] √ π cosh (π ), (29)
where

k 1 (a, σ * t ) = √ 2 π a * a σ * t (x) sin ( ln (x -a)) √ x -a dx, k 2 (a, σ * t ) = √ 2 π a * a σ * t (x) cos ( ln (x -a)) √ x -a dx, (30)
and the parameter accounts for the different properties of the two bonded solids, in case one of them being rigid, it equals to

= - 1 2π ln (3 -4ν). (31)
In Fig. 14(d) we plot the approximation of the shear traction distribution in the interface during trap opening, discussed above. A sound similarity is found between numerical results and analytical formulae provided by the LEFM. Therefore, we have shown that during the process of trap opening due to increasing pressure in the fluid with friction taken into account, the tangential tractions near the contact edges are elevated up to the limit provided by the Coulomb friction law. Consequently, even if the majority of the interface remains in stick state, local slip zones emerge at the boundaries of contact zones. It is important to account for such an elevated shear stress near edges of contact zones surrounding trapped fluid, which appears under normal loading, in the analysis of damage evolution and crack onset under monotonic and cycling loading, including fretting fatigue [START_REF] Hills | Mechanics of fretting fatigue[END_REF][START_REF] Proudhon | A fretting crack initiation prediction taking into account the surface roughness and the crack nucleation process volume[END_REF].

Conclusions

In this work we solved the problem of mechanical contact between a deformable body with a wavy surface and a rigid flat, taking into account pressurized fluid trapped in the interface. A mathematical framework for this coupled problem for both incompressible and compressible fluids was formulated. In the latter case, either constant or pressure-dependent fluid bulk-moduli were considered; all models were implemented in the finite element framework using a monolithic approach.

The proposed framework accounts for a finite slope of the roughness profile, while in previous investigations using classical boundary element method (which accounts only for vertical displacements) and existing analytical solutions only infinitesimal slopes were considered. We show that in the considered coupled problem, a reduction of the contact area can occur due to elastic flattening of asperities by fluid pressure. Thus the reduction of the global coefficient of friction is caused not only by the external load repartition between the solid contact and the pressurized fluid, but also by the contact area reduction.

The reduction of the contact area takes place if the fluid bulk-modulus is higher than that of the solid. In case of incompressible fluid this criterion is satisfied and the process of trap opening is observed. However, this case is non-physical, since real lubricating fluids in the unpressurized state have much lower bulk modulus than solids. A more relevant case is a compressible fluid with linear dependence of bulk modulus on pressure, which ensures a non-monotonic variation of the contact area, and thus of the global coefficient of friction, leading to reduction of the both for sufficiently large pressures.

Among other applications, the obtained results are relevant for the mechanical behaviour of multi-cracked materials such as rocks. We showed that due to the presence of pressurized fluid in the interface, the frictional limit does not increase linearly with increasing external load, but reaches its maximum and decreases down to zero. This behaviour is similar to pressure-dependent plasticity models with a cap (e.g. Drucker-Prager cap model), which corresponds to the decay of the von Mises yield stress with the increasing pressure.

In addition to elasticity, we considered physically more relevant elasto-plastic materials in combination with realistic fluids. In this case, the contact pressure is bounded, while the fluid can bear arbitrary pressure, consequently under certain external pressure fluid permeates in the contact zones abruptly.

When interfacial friction is considered in the coupled problem, previously unreported quasi-singularities appear in shear stresses near edges of contact patches during fluid-trap opening under normal loading. We showed that these singularities can be analytically estimated using the analogy between trap opening and crack propagation in the interface between two bonded dissimilar solids. It is important to account for such an elevated shear stress, caused by the trapped fluid, in the analysis of damage evolution and crack onset under monotonic and cycling loading, including fretting fatigue.

The problem of trapped fluid is relevant for metal forming (drawing and rolling), where a lubricant is present in the interface and involved loads are high. It is also relevant in poromechanics, especially in cracked media filled with fluid and subjected to complex stress states with high hydrostatic component, which can ensure contact between surfaces of internal cracks. Finally, at the microscopic scale, where the surface roughness plays a crucial role, the trapped fluid provides additional load-bearing capacity, and thus reduces the macroscopic static friction. Under increasing load, the trapped fluid is squeezed out of its trap thus resulting in even smaller global coefficient of friction. anonymous reviewer are kindly acknowledged. According to the integrated Flamant's solution, any uniformly distributed pressure on the surface will result in a uniform vertical displacement. However, it is true only if the surface is flat. For a wavy surface, under the action of a uniform pressure, the crest's displacement is bigger than the displacement of the trough. The uniform pressure distribution on the surface is given by -p 0 n, where n is the outer normal to the surface (see Fig. A.15). We consider the vertical p and horizontal q components of the normal pressure (each one of them contributes to the distortion of the profile):

q(x) = -2πp 0 ∆ λ sin 2πx λ + O ∆ 3 λ 3 and p(x) = -p 0 + O ∆ 2 λ 2
in case of small slope. Therefore, keeping the small values of order ∆/λ, we may calculate the vertical displacement caused by the horizontal component q using the integrated Flamant's solution:

u q y (x) = - (1 -2ν)(1 + ν) 2E    x -b
q(s) ds -b x q(s) ds

   + C,
where b → ∞ and C is an arbitrary constant. Substituting q(s) and calculating integrals, we obtain:

u q y (x) = - (1 -2ν)(1 + ν)p 0 E ∆ cos 2πx λ + C. (A.1)
We are also convinced that the uniformly distributed vertical traction p produces exactly the same vertical displacement:

u p y (x) = - (1 -2ν)(1 + ν)p 0 E ∆ cos 2πx λ + D, (A.2)
where D is another arbitrary constant. However, this result does not follow from Flamant's solution, which as already mentioned would predict a uniform displacement. This result was guessed and confirmed with a very high accuracy (fractions of percent) by finite-element simulations for different fractions ∆/λ and Poisson's ratios. Of course, the simplest analogy would be the Winkler's foundation with springs whose lengths follow the distribution l(x) = L + ∆ cos(2πx/λ), which would mimic the shape of the wavy surface. However, it is unavailing to obtain the proportionality factor of form (1 -2ν)(1 + ν)/E. At this stage we are content with numerical proof only, which consisted in applying separately horizontal and normal components of the pressure over a single period with periodic boundary conditions and comparing the numerical results with the equations (A.2) and (A.1). The total displacement field reads:

u y (x) = u q y (x) + u p y (x) = - 2(1 -2ν)(1 + ν)p 0 E ∆ cos 2πx λ + C.
We may define C so that u y (0) = y 0 , where y 0 is the current vertical position of the crest, thus giving:

u y (x) = y 0 - 2(1 -2ν)(1 + ν)p 0 E ∆ cos 2πx λ -1 .
We calculate the volume between the distorted wavy surface and the plane y = y 0 :

V g = Γ (y 0 -(Y (x) + u y (x)))dΓ,
where Y (x) = ∆(cos 2πx λ -1) is the initial vertical coordinate. Note that horizontal displacements do not contribute significantly to the volume change. Upon integration, noting that the initial gap volume V g0 = ∆λ, we obtain:

V g /V g0 = 1 - 2(1 -2ν)(1 + ν)p 0 E ,
which gives the relative change of the volume between the deformed wavy surface and corresponding plane under action of the uniform hydrostatic pressure p 0 . Note that it does not depend on ∆/λ, however, this result was obtained under assumption of ∆/λ 1, therefore it corresponds to the limiting case of small, but finite slope of the profile. In order to prove validity of this extension, we will consider a transformation of the Lagrangian for the coupled system (17). We will start with substituting the formula for the gap volume (15) into ( 17) and obtain:

Γ f Γ * Γ c Γ c Γ *
L(u, λ c , λ f ) = Π(u) + Γc λ c g(u) dΓ c -λ f    Γ f g(u) d Γ f -V f    . (B.8)
Let us break the integral over the active contact zone Γ c into two integrals over surfaces Γ * and Γ c \ Γ * and consider the following representation of contact Lagrange multiplier λ c on the surface Γ * : λ c = λ * c -λ f , where λ * c ≤ 0. Note that this representation is valid only if |λ c | ≥ |λ f |, which is the case for the problem under study (except for the elasto-plastic case): due to the considered regular wavy profile of the surface and gradual monotonic increase of the external pressure, the contact pressure must be higher in the contact patches (i.e. in the active contact zone), than in the trapped fluid zone, because otherwise the contact would not be active. Therefore, we may write

L(u, λ c , λ f ) =Π(u) + Γc\Γ * λ c g(u) dΓ c + Γ * (λ * c -λ f ) g(u) dΓ c -λ f    Γ f g(u) d Γ f -V f    =Π(u) + Γc\Γ * λ c g(u) dΓ c + Γ * λ * c g(u) dΓ c -λ f    Γ f ∪Γ * g(u) d Γ f -V f    . (B.9)
The last formula in (B.9) shows that for the coupled problem:

• the trapped fluid zone Γ f can be extended on a part or on the whole active contact zone Γ c without loss of generality: Γ * ⊆ Γ c ;

• if the trapped fluid is in the active state, the value of λ * c -λ f is equivalent to the normal stress component on Γ * .

  Figure 2: A sketch of the problem under study.

Figure 3 :

 3 Figure 3: (a) Reference configuration: X ∈ Ω 0 . (b) Actual configuration: x ∈ Ω, p 0 is the external pressure.

Figure 4 :

 4 Figure 4: (a) Trapped fluid is in inactive state. (b) Trapped fluid is in active state, p f is the fluid pressure.

Figure 5 :

 5 Figure 5: FEM mesh

  The local Coulomb's coefficient of friction determines the following inequality: |σ t | ≤ µ loc |σ n |, where σ t and σ n are the tangential and normal components of the traction vector in the contact interface, respectively.

FEMFigure 6 :

 6 Figure6: (a) The evolution of the real contact area in the vicinity of the "activation" point of the incompressible fluid with respect to the external pressure normalized by E * for three profiles with different slopes ∆/λ and three cases with different ratios of the fluid volume to the initial gap volume V f /V g0 . (b) The evolution of the real contact area until the complete opening of the trap for the case V f /V g0 = 0.9 shown for different slopes ∆/λ.

Figure 7

 7 Figure 7: (a) Sketch of the auxiliary problem: deformation of the wavy surface under uniform hydrostatic pressure. (b)Evolution of the ratio V f /V g0 (V f is the volume between the deformed surface and a horizontal plane y = y 0 , where y 0 is the current position of the crest, and V g0 is the initial volume of the gap) with the increasing hydrostatic pressure p f for several profiles with different slope ∆/λ.

Figure 8 :

 8 Figure 8: Stress and strain components in the bulk of the deformable solid during the process of trap opening due to the increasing pressure in the fluid. Top to bottom: vertical stress component σyy, von Mises stress σ vM , hydrostatic stress p, horizontal strain component εxx and the vertical one εyy. Three loading steps are considered, corresponding to, left to right: maximal contact area (activation of the fluid), half of the contact opened, contact area is zero (trap is opened). The considered elastic material is typical aluminium (E = 70 GPa, ν = 0.33), the fluid is assumed incompressible.

Figure 9 :

 9 Figure 9: (a) Real contact area evolution during opening of the contact, caused by pressurized incompressible trapped fluid (with respect to the external pressure normalized by E * ) (b) Contact normal force evolution during opening of the contact (with respect to the external pressure, normalized by E * ). (c) Evolution of the ratio between global and local coefficients of friction, and (d) a zoom of this evolution for V f /V g0 = 0.9, where, in addition, the results of frictional simulations are plotted (crosses), as well as analytical approximations given by (27) (dashed curves).

  , p 0 / p * Normalized real contact area, A/A 0

Figure 10 :

 10 Figure10: (a) Evolution of the ratio of the real contact area to the apparent one under increasing external pressure p 0 : comparison of numerical (solid curves with markers) and analytical (dashed curves) results for different values of the fluid modulus of compressibility, normalized by the bulk modulus of the solid K f /Ks; ∆/λ = 0.01, V f 0 /V g0 = 0.7. (b) Distribution of the normal pressure near the contact patch under the increasing external load p 0 . Solid lines are the results of the numerical simulation and dashed lines correspond to the analytical solution under the same external pressure, ∆/λ = 0.01, V f 0 /V g0 = 0.9, K f /Ks = 6 • 10 4 .

Figure 11 :

 11 Figure11: Evolution of (a) the ratio of real contact area to the apparent one, (b) the ratio between global and local coefficients of friction under increasing external pressure for two elastic solids representing steel and aluminium, and non-linearly compressible fluids representing water, glycerine and oil. The dashed curves correspond to the analytical solution given by (12). Evolution of (c) the ratio of real contact area to the apparent one A/A 0 , and of (d) the global to local coefficients of friction under increasing external pressure in the case of elastic-perfectly plastic solid and incompressible fluid. Note that in the initial configuration the fluid does not occupy the entire gap: V f /V g0 = 0.9. Dashed curves are presented for comparison with the cases of purely elastic solids, discussed in Sec. 5.1. Vertical dash-dotted line indicates the hardness taken to be H = 3σ Y .

Figure 12 :

 12 Figure 12: Accumulated plastic strain near the contact interface is shown at three different external loads for the incompressible fluid and for V f /V g0 = 0.9 and ∆/λ = 0.01, from left to right: (1) the step corresponding to activation of the fluid; (2) the plastic zone appears in the valley between the asperities; (3) opening of the trap (contact area is zero; as seen from the plastic field, at this moment the entire solid is plastified).

Figure 13 :Figure 14 :

 1314 Figure 13: The behaviour of the system considering elasto-plastic material and non-linearly compressible fluid: (a) evolution of the ratio of real contact area to the apparent one under increasing external pressure; (b) the same as (a), but the results are shown in range 0 ≤ p 0 ≤ 0.025E * ; (c) evolution of the ratio between global and local coefficients of friction; (d) the same as (c), but the results are shown in range 0 ≤ p 0 ≤ 0.025E * . Dashed curves are presented for comparison with the cases of purely elastic solids. Vertical dash-dotted line indicates the hardness H = 3σ Y .

  Figure A.15: Distortion of a periodic wavy surface under a uniform hydrostatic pressure.

Figure B. 17 :

 17 Figure B.17: Extension of the trapped fluid zone Γ f on the active contact zone Γc.

  map:81.0000 time:0.81 min:-10282.3 max:-10077.0

																	highlights the contact zone
					rigid flat											rigid flat										rigid flat
	-7.2e+02 -0.72 -6.6e+02	-6e+02	-5.4e+02	-4.7e+02	-4.1e+02 sig22 map:9.00000 time:0.09 -3.5e+02 min:-720.309 max:15.9072 -2.9e+02 σ yy (GPa)	-2.3e+02	-1.7e+02	-1.1e+02	-45	16 0	-7e+03 -7	-7e+03	-6.9e+03	-6.9e+03	-6.9e+03	-6.8e+03 sig22 map:67.0000 time:0.67 -6.8e+03 min:-7006.45 max:-6564.98 -6.7e+03 σ yy (GPa)	-6.7e+03	-6.7e+03	-6.6e+03	-6.6e+03 -6.6 -6.6e+03	-10.3			σ yy (GPa)	-10.1
						rigid flat											rigid flat										rigid flat
	0 -7.2e+02	-6.6e+02	-6e+02	-5.4e+02	-4.7e+02	-4.1e+02	-3.5e+02	-2.9e+02	-2.3e+02	-1.7e+02	-1.1e+02	0.41 3.1 16 -45 -7e+03	-7e+03	-6.9e+03	-6.9e+03	-6.9e+03		-6.8e+03	-6.8e+03	-6.7e+03	-6.7e+03	-6.7e+03	-6.6e+03	-6.6e+03	3.6 4.8 -1e+04 -1e+04 -6.6e+03	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	5.6 -1e+04
	0 -7.2e+02	-6.6e+02	-6e+02	-5.4e+02	-4.7e+02	-4.1e+02	-3.5e+02	-2.9e+02	-2.3e+02	-1.7e+02	-1.1e+02	0.47 3.2 16 -45													3.7 4.9 -1e+04 -1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	-1e+04	5.3 -1e+04
						rigid flat											rigid flat										rigid flat
	-0.7% -7.2e+02 -6.6e+02 -6e+02	-5.4e+02	-4.7e+02	-4.1e+02 sig22 map:9.00000 time:0.09 -3.5e+02 min:-720.309 max:15.9072 -2.9e+02 ε xx		-2.3e+02	-1.7e+02	0.7% 16 -1.1e+02 -45	-7e+03 -0.7% -7e+03 -6.9e+03	-6.9e+03	-6.9e+03	-6.8e+03 sig22 map:67.0000 time:0.67 -6.8e+03 min:-7006.45 max:-6564.98 -6.7e+03 ε xx	-6.7e+03	-6.7e+03	-6.6e+03	-6.6e+03 0.7% -6.6e+03	-1e+04 -0.7% -1e+04 -1e+04	-1e+04	-1e+04	-1e+04 sig22 map:81.0000 time:0.81 -1e+04 min:-10282.3 max:-10077.0 -1e+04 ε xx	-1e+04	-1e+04	-1e+04	-1e+04 0.7% -1e+04

sig22 map:81.0000 time:0.81 min:-10282.3 max:-10077.0 σ vM (GPa) sig22 map:9.00000 time:0.09 min:-720.309 max:15.9072 σ vM (GPa) sig22 map:67.0000 time:0.67 min:-7006.45 max:-6564.98 σ vM (GPa) sig22 map:81.0000 time:0.81 min:-10282.3 max:-10077.0 p (GPa) sig22 map:9.00000 time:0.09 min:-720.309 max:15.9072 p (GPa)

As was mentioned earlier, the infinitesimal-slope assumption implies here that (i) the wavy surface behaves as a flat one and that Flamant's solution holds for every surface point, and (ii) that the horizontal component of the fluid pressure is negligible.

contact zone, i.e. consider fluid pressure on surface Γ f ∪ Γ * , where Γ * ⊆ Γ c , see Fig.B.17. The only change we have to take into account is that on Γ * the contact normal pressure will not be equal to the Lagrange multiplier λ c corresponding to contact, but to the difference between the latter and the value of the trapped fluid Lagrange multiplier: σ n = λ c -λ f on Γ * . Note that λ c , which is equivalent to the normal traction, is negative, while λ f represent fluid pressure, which is positive by definition.

We demonstrate it for the case of the Lagrange multipliers method, however, it may be also generalized to the penalty method.

Acknowledgements

The authors acknowledge the financial support of Safran Tech and MINES ParisTech (Thèse-Open) and are grateful to Julien Vignollet for his helpful suggestions and kind support. Enlightening comments of an

Appendix B. Trapped fluid element

In order to implement in a finite element code our approach for modelling the trapped fluid, alongside with structural and contact elements, we used a special trapped-fluid element containing all segments of the trapped fluid zone Γ f , see Fig. B.16 (in the FEM literature this element is also known as a hydrostatic fluid element).

In the finite element framework the area of the gap (15) can be calculated by the following formula:

where the summation is performed over all segments of the surface Γ f , X y i and u y i are the vertical coordinate in the reference configuration and the vertical displacement of the i-th node of the corresponding segment, respectively, N i (ξ) is the shape function, associated with the i-th node; n y seg (ξ) = ν • n seg (ξ), where ν is the normal to the rigid plane, and n seg (ξ) is the normal to the segment, J(ξ) is the Jacobian and ξ ∈ [-1, 1] is the convective coordinate in the parent space. Note that summation over the repeating indices is assumed.

Therefore we can consider the gap volume (B.1) as a function of the displacement vector u = [u x 1 , u y 1 , . . . u x N , u y N ] T , consisting of the displacement components of N nodes on the surface Γ f , (u x i , u y i ) are horizontal and vertical components of the displacement vector of the i-th node, respectively.

Appendix B.1. Lagrange multiplier formulation

The problem of finding a stationary point of the Lagrangian ( 17) is non-linear, and to solve it numerically we use the classical Newton-Raphson method, which requires calculation of the residual vector of the trappedfluid element

and the corresponding tangent matrix

Note that the vector of degrees of freedom (DOFs) for the trapped fluid element has the form [u x 1 , u y 1 , . . . u x N , u y N , λ f ] T , and includes 2N displacement DOFs of N surface nodes and also an additional single Lagrange multiplier λ f (which represents the fluid pressure).

On every global iteration of the Newton-Raphson method we perform the following active set strategy, see [START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF] for the trapped fluid element:

1. Get from previous iteration current value of the displacement vector u and the Lagrange multiplier λ f , calculate the gap volume V g (u); 2. If the gap V g (u) > V f or λ f < 0, the fluid is considered to be in inactive state and the trapped fluid element is excluded from consequent calculations; 3. Otherwise, i.e. if V g (u) ≤ V f and λ f ≥ 0, the fluid is in active state, then the residual vector and the tangent matrix are calculated by formulas (B.2) and (B.3).

Appendix B.2. Penalty formulation

For the numerical simulations of the trapped compressible fluid we may consider the same finite element, as was described previously for the case of Lagrange multipliers, with one difference: no extra degrees of freedom are involved, and the vector of DOFs has the form [u x 1 , u y 1 , . . . u x N , u y N ] T . A similar active set strategy can be applied to this element, as it was described for the case of Lagrange multiplier element.

The residual vector and the tangent matrix for the linear penalty trapped fluid element have the form:

where ⊗ is a tensor product, see [START_REF] Yastrebov | Numerical methods in contact mechanics[END_REF]. In the case of a non-linear penalty element: Initially we supposed that active contact zone Γ c and trapped fluid zone Γ f are complementary subsets of Γ -the whole surface of the deformable body (in the interface): Γ f ∩ Γ c is a set of measure zero, and Γ f ∪ Γ c = Γ. In accordance with the numerical procedures for solving the coupled problem proposed in the previous subsections, the trapped fluid zone Γ f , and, consequently, the number of DOFs of the trapped fluid element must be updated on every iteration of the Newton-Raphson method, which increases the computation time. Below we show that in order to simplify numerical calculations we may omit this zone splitting without loss of generality and accuracy 2 . We make an extension of trapped fluid zone on the active