Computational framework for monolithic coupling for thin fluid flow in contact interfaces

by   Andrei G. Shvarts, et al.

We develop a computational framework for simulating thin fluid flow in narrow interfaces between contacting solids, which is relevant for a range engineering, biological and geophysical applications. The treatment of this problem requires coupling between fluid and solid mechanics equations, further complicated by contact constraints and potentially complex geometrical features of contacting surfaces. We develop a monolithic finite-element framework for handling contact, thin incompressible viscous flow and fluid-induced tractions on the surface of the solid, suitable for both one- and two-way coupling approaches. Additionally, we consider fluid entrapment in "pools" delimited by contact patches and its pressurisation following a non-linear compressible constitutive law. Image analysis algorithms are adopted to identify the local status of each interface element (i.e. distinguish between contact, fluid flow and trapped fluid zones) within the Newton convergence loop. First, an application of the proposed framework for a problem with a model geometry is given, and the robustness is demonstrated by the DOF-wise and status-wise convergence. The full capability of the developed two-way coupling framework is demonstrated on a problem of a fluid flow in a contact interface between a solid with representative rough surface and a rigid flat. The evolution of the contact pressure, fluid flow pattern and the morphology of trapped fluid zones under increasing external load until the complete sealing of the interface is displayed. Finally, effective properties of flat-on-flat rough contact interfaces such as transmissivity and real contact area growth are calculated using the developed framework, showing qualitatively new results compared to the one-way coupling approximation.



There are no comments yet.


page 10

page 13

page 15

page 20

page 21

page 23

page 26

page 27


A consistent approach for fluid-structure-contact interaction based on a porous flow model for rough surface contact

Simulation approaches for fluid-structure-contact interaction, especiall...

Frontiers in Mortar Methods for Isogeometric Analysis

Complex geometries as common in industrial applications consist of multi...

3D-2D Stokes-Darcy coupling for the modelling of seepage with an application to fluid-structure interaction with contact

In this note we introduce a mixed dimensional Stokes-Darcy coupling wher...

Fabrication of Embedded Microvalve on PMMA Microfluidic Devices through Surface Functionalization

The integration of a PDMS membrane within orthogonally placed PMMA micro...

The unstable temporal development of axi-symmetric jets of incompressible fluid

We study the shear-driven instability, of Kelvin-Helmholtz type, that fo...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The problem of thin fluid flow in narrow interfaces between contacting or slightly separated deformable solids appears in different contexts: in engineering and biological applications, as well as in geophysical sciences. Rigorous handling of such problems requires solution of a strongly nonlinear contact problem, which is further complicated by a multi-field coupling of essentially interrelated fluid and solid mechanics. The free volume111By free volume, here, we mean the void separating the contacting surfaces. between contacting surfaces depends on their initial geometry, which can be rather complex: they may have deterministic features (e.g. turned (Pérez-Ràfols et al., 2016), patterned surfaces (Sahlin et al., 2010b; Prodanov et al., 2013)) and, at a certain magnification, may be considered as randomly rough (Nayak, 1971) down to atomistic scale (Krim and Palasantzas, 1995; Ponson et al., 2006; Whitehouse, 2010; Thomas, 1999).

Numerous applications of the problem of thin fluid flow in contact interfaces include sealing engineering (Müller and Nau, 1998), lubrication in elasto-hydrodynamic and mixed regimes (Stupkiewicz and Marciniszyn, 2009; Sahlin et al., 2010a) and functioning of human joints (Caligaris and Ateshian, 2008). Such an interaction between fluids and solids in contact is also relevant for hydraulic fracturing (Bažant et al., 2014), extraction of shale gas and oil from rocks (Hubbert and Willis, 1972), and at larger scales in landslides (Viesca and Rice, 2012), slip in pressurized faults (Garagash and Germanovich, 2012) and bazal sliding of glaciers (Fischer and Clarke, 1997).

The problem under discussion belongs to a vast domain of fluid-structure interaction (FSI) problems, which involve deformation and/or motion of the solids interacting with the internal and/or external fluid(s). These problems are of very wide range, spanning from deformation of airplane wings and rotor blades subjected to the sub- or supersonic air flow (Farhat et al., 2003; Bazilevs et al., 2011) to modelling of the blood flow (Bazilevs et al., 2006; Gerbeau and Vidrascu, 2003) and heart valves (Dos Santos et al., 2008; Astorino et al., 2009; Kamensky et al., 2015), scaling up to suspended bridge instabilities under wind load (Païdoussis et al., 2010), ship stability (Wackers et al., 2011) or huge iceberg’s capsize in water (Sergeant et al., 2018). All these problems correspond to different space and time scales, operating conditions and other requirements, therefore a unified FSI approach fit for all cases does not exist, and different techniques have been developed for particular problems.

Many problems of the fluid-structure interaction, such as aeroelasticity and hemodynamics, correspond to the case of the high-Reynolds-number flow. Therefore, different mesh density, and often different time stepping, are required for the solid and fluid domains. Furthermore, the fluid domain evolves due to motion and deformation of solids. A number of methods have been used to overcome the associated computational complexity, such as arbitrary Lagrangian-Eulerian method (Donea et al., 1982; Takashi and Hughes, 1992; Wick, 2014), fictitious domain method (Baaijens, 2001; De Hart et al., 2003; Dos Santos et al., 2008), immersed boundary method (Mittal and Iaccarino, 2005; Kadapa et al., 2017) and extended finite element method (Mayer et al., 2010; Gerstenberger and Wall, 2008). On the contrary, fluid flow in contacting or slightly separated interfaces is usually of low Reynolds number and, moreover, the thickness of the fluid film is usually much smaller than other dimensions of the solid. In this case general Navier-Stokes equations could be readily simplified down to the Reynolds equation for the viscous flow (Hamrock et al., 2004). This simplification permits to use compatible meshes for the fluid and the solid domains and, under assumption of constant pressure through the film thickness, define the Reynolds equation on the so-called lubrication surface, so that specific methods discussed above are not required, see (Stupkiewicz and Marciniszyn, 2009; Stupkiewicz et al., 2016).

From the point of view of the underlying physical processes, different FSI strategies could be divided into one-way and two-way coupling approaches. In the context of thin fluid flow through contact interfaces, the former implies that the solution of the solid mechanics problem defines the distribution of the free volume in the interface which can be occupied by the fluid flow, however the fluid pressure does not affect the deformation of the solids, i.e. the fluid problem is solved assuming rigid walls of solids. In the two-way coupling this approximation is dropped, and the fluid-induced traction acting on the surface of the deformable solid is taken into account.

In elastohydrodynamic lubrication regime, as well as for non-contact seals, the two-way coupling is often used, see (Stupkiewicz and Marciniszyn, 2009; Stupkiewicz et al., 2016; Yang and Laursen, 2009). However, for the important case of contact seals, or, more generally, if contact is present in the interface, the one-way coupling is rather utilized, see (Dapp et al., 2012; Dapp and Müser, 2016; Pérez-Ràfols et al., 2016)

. It is widely assumed in this context that the deformation of the solids happens mainly due to the contact interaction, and the fluid pressure effect on the solid is negligible, since the contact pressure at surface’s asperities is considerably higher than the physically relevant fluid pressure. However, to the best of authors knowledge, a quantitative range of validity of one-way coupling for problems involving thin fluid flow in contact interfaces, depending on the surface geometry, material properties and the fluid pressure, has not been determined yet. The lack of such quantitative analysis is probably caused by insufficiency of existing numerical methods for comparison of one- and two-way coupling approaches for these problems, which is the main motivation of the current study.

It is important to note, that the application of the above mentioned rather general FSI approaches to problems involving solid-solid contact provide possibility for two-way coupling, see e.g. (Mayer et al., 2010; Kamensky et al., 2015). However, these methods become inefficient if a very fine discretization is used to accurately represent surface roughness in a numerical simulation. This computational complexity can be overcome by considering surface roughness in an averaged sense at the macroscopic level, see (Patir and Cheng, 1978; Pérez-Ràfols et al., 2016; Zaouter et al., 2018; Waseem et al., 2017), or as was done recently using a porous flow model (Ager et al., 2018). Nonetheless, if one is interested in the effect of rough or deterministic surface features on the solution of the coupled problem, e.g. distribution of contact stresses and fluid flow patters, this approximation is not appropriate. Therefore, the main purpose of the current study is to develop a computational framework suitable for both one- and two-way coupling approaches and applicable for a discretization which reflects essential features of the surface geometry.

An important example of the effect of the surface roughness on the coupled problem, relevant only for the two-way coupling, is the phenomenon of fluid entrapment. Indeed, studying the evolution of the morphology of the contact interface under increasing normal load, one may observe how non-simply connected contact patches appear, see (Yastrebov et al., 2015; Lorenz and Persson, 2010). At the same time, the fluid present in the interface can be trapped inside “pools” or “valleys” bounded by these patches and subsequently become pressurized and provide additional load-bearing capacity. The behaviour of trapped fluid accounts for a significant reduction of friction in tire-road contact (Scaraggi and Persson, 2012), cold metal forming (Azushima and Kudo, 1995) and functioning of human joints (Soltz et al., 2003; Chan et al., 2011). However, the effect of the fluid entrapment on transmissivity of contact interfaces, which is an import characteristic for sealing engineering, has not been thoroughly investigated yet.

We have recently studied the behaviour of hydrostatically pressurized fluid trapped in a contact interface without considering the fluid flow (Shvarts and Yastrebov, 2018b)

. We used the augmented Lagrangian method for the contact constraints, and the classic Lagrange multiplier method permitted us to take into account the additional constraint in case of incompressible fluid, while the penalty method was applied if the fluid was assumed compressible. We introduced a novel trapped fluid super-element based on all out-of-contact segments (faces in 3D), surrounded by a contact patch, so that displacements of all nodes associated with this element were included into its degrees of freedom (DOF) vector. Moreover, we proposed a technique of extension of the trapped fluid element on the surrounding contact zone by superposing locally the contact and fluid pressure fields, which simplified considerably computations in certain cases (see Appendix B in 

(Shvarts and Yastrebov, 2018b)). Therefore, an additional objective of the current study is to incorporate the formulation of the trapped fluid element into the computational framework, permitting to quantify the effect of trapped fluid zones on the coupled problem.

From the implementation point of view two distinct approaches for any FSI problem exist: partitioned and monolithic. The former is based on two different solvers for the fluid and solid sub-problems, and in order to take into account the coupling, one- or two-way data exchange between them must be established. Furthermore, a certain iterative process is required to obtain the convergence. The utilization of the partitioned approach benefits from modularity, since different solvers tailored for the sub-problems could be used (Küttler and Wall, 2008; Matthies and Steindorf, 2003), however, convergence and stability of such scheme could raise issues, and special techniques may become necessary, see, for example (Heil, 1998). On the contrary, under the monolithic approach all equations which govern sub-problems and the interaction between them are rendered into one system, and upon its solution DOF values corresponding to both sub-problems are obtained simultaneously (Hübner et al., 2004; Michler et al., 2004; Heil, 2004; Verdugo and Wall, 2016). The data exchange in this case is not needed, the stability and convergence are obtained easier, however, solution of the vast system of algebraic equations is necessary. Nevertheless, for the problem under discussion this last issue is not relevant, since the number of unknowns in the interface is considerably smaller than in the bulk of the solid, if, for example, finite-element discretization is used.

Elasto-hydrodynamic lubrication problems are often solved under the monolithic approach (Stupkiewicz et al., 2016; Stupkiewicz and Marciniszyn, 2009; Yang and Laursen, 2009), whilst for the contact sealing problems the partitioned approach is generally preferred (Pérez-Ràfols et al., 2016). Moreover, as was already mentioned above, the problem is often solved under the one-way coupling approach, using the assumption of the infinitesimal slopes of the surface profile and the small deformation formulation. Boundary element method (Pérez-Ràfols et al., 2016) and Green’s function molecular dynamics (Dapp et al., 2012) are frequently used for the mechanical contact problem and the Reynolds equation is often solved by the finite-differences method.

In this paper we develop a computational framework aimed first at resolution of the two-way coupling of the mechanical contact and fluid flow sub-problems, which could require rather frequent and considerable data exchange in case of the partitioned approach. The second objective is to take into account the effect of trapped fluid pools appearing in the interface. If the number of these zones becomes large, then a resolution under the partitioned approach becomes even more complicated, since the history tracking of trapped zones is needed. The monolithic approach appears beneficial for our purposes and, therefore, it was applied throughout this study. We use the finite-element method in order to make possible application of the proposed framework for different surface geometries (e.g. with finite slopes of the profile), under larger deformation formulation and with different material models of the solid (e.g. elasto-plastic, viscoelastic, etc). The Newton-Raphson method is applied to solve the essentially nonlinear problem, and, furthermore, on each iteration we perform identification of the local status in the interface to distinguish between contact zones and fluid flow zones and also keep track of formation and evolution of trapped fluid pools.

Note that the methodology presented here was already applied to a two-way coupling problem of a pressure driven fluid flow in contact interface between an elastic solid with an extruded wavy surface and a rigid flat (Shvarts and Yastrebov, 2018a). In the same paper we derived an approximate analytical formula based on the Westergaard-Kuznetsov solution (Westergaard, 1939; Kuznetsov, 1985) and a one-dimensional formulation of Reynolds equation, which describes both the solid deformation and the fluid pressure distribution in case of two-way coupling. The comparison of this approximation with a numerical solution based on the presented framework showed good agreement between the two in the interval of loads within which the analytical one is applicable. Furthermore, the numerical method permitted us to investigate the problem with parameters beyond the range of validity of the analytical solution up to full sealing of the interface. However, since only regular wavy profiles were studied, the problem did not not include trapped fluid zones, therefore a complete testing of the proposed method was not possible. Moreover, the methodology was only used, and was not discussed in detail in the above cited work.

The present paper is organized as follows. In Section 2 we formulate the coupled problem by outlining equations governing each sub-part of it, while Section 3 introduces the variational statement for this problem under both one- and two-way coupling approaches. In Section 4 we propose the monolithic framework and provide ready-to-implement finite-element formulations of the tangent matrix and the residual vector of the coupled problem. Moreover, we discuss the algorithm used to determine the local status in the interface and keep track of formation and evolution of multiple trapped fluid zones. Section 5 is devoted to examples showing capabilities of the proposed framework and some relevant discussions. Finally, Section 6 provides a short conclusion.

2 Problem statement

We consider the problem of a thin fluid flow in contact interface between a solid body with an arbitrary surface geometry, given, for concreteness, by a function , and a rigid flat described by a plane . This set-up was chosen to simplify the formulation of the contact problem and to concentrate the discussion on handling of the fluid/solid coupling. Note that the unilateral contact set-up is equivalent to the case of contact between two deformable bodies under linear isotropic elasticity and infinitesimal strain assumptions, see (Barber, 2013). Furthermore, the problem statement with one of the contacting solids considered as rigid is relevant when the elasticity modulus of this solid is considerably higher than of the other, which is the case, for example, in rubber sealing applications (Persson and Yang, 2008) and in tire-road contact (Scaraggi and Persson, 2012). At the same time, the numerical framework developed in this paper for the case of unilateral contact can be extended afterwards to the problem of contact between two deformable solids with arbitrary surface geometries.

Let us denote by the deformable solid and by the part of its surface which represents the potential contact zone, i.e. defines the extent of the contact interface. The resolution of the coupled problem requires subdivision of the surface into following parts, according to the local status222By the local status in this context we mean the location of each point of the solid’s surface in the active contact zone, fluid-structure interface, or one of the trapped fluid zones. of the interface, see Fig. 1:


where is the active contact zone where normal contact tractions are non-zero, is the part of the solid’s surface which interacts with the flowing fluid and where the surface tractions are equal to the corresponding tractions in the fluid (so-called fluid-structure interface), are trapped fluid zones, i.e. parts of the surface which are out of contact, but completely delimited by non-simply connected contact patches. Furthermore, we term by the projection of on the plane , which serves as the lubrication surface where the Reynolds equation for the fluid flow will be defined.

By definition, and . Note also, that , i.e. even though all trapped pockets contain the same fluid, as the one present in the fluid-flow domain, the behaviour of the trapped fluid is to be considered separately. We assume here that in the initial configuration and , so that the fluid flow zone occupies the whole interface. This assumption makes impossible appearance of non-contact zones which do not belong to or one of the trapped fluid zones in any deformed configuration. However, consideration of such a case (e.g. air bubbles entrapment) could be included into the framework.

Figure 1: Sketch of the problem under study: (a) contact between a solid and a rigid flat with fluid present in the interface; (b) view of the contact interface. Notations: is the active contact zone, is the fluid-structure interface, is the lubrication surface where the Reynolds equation is solved ( is the projection of on the rigid flat), is a trapped fluid zone.

2.1 Solid mechanics problem with unilateral contact

The deformation of the solid (in absence of the fluid) is governed by the balance of momentum equation complemented by the contact and boundary conditions:

in (2a)
on (2b)
on (2c)
on (2d)

where is the displacement field ( and are the coordinates in the deformed and the initial configurations, respectively),

is the Cauchy stress tensor; (

2b) are the Hertz-Signorini-Moreau conditions of the non-adhesive frictionless unilateral contact (Wriggers, 2006; Yastrebov, 2013), is the normal traction on the solid’s surface (to which is the outer normal). We denote by the normal gap function, i.e. the signed distance between the surface of the solid and the rigid flat: in case of separation, in contact and refers to penetration, which is not admissible:


where is the initial gap and is the normal to the rigid flat. Finally, (2c) are the Dirichlet boundary conditions with a prescribed displacement and (2d) are the Neumann boundary conditions with a prescribed surface traction , defined on and , respectively. Note that here we will consider isotropic linearly elastic solid, i.e. the Cauchy stress tensor is related to the infinitesimal strain tensor by the Hooke’s law: with the Lamé constants (elastic moduli) and , is the identity tensor. However, the proposed coupled framework concerns processes occurring in the contact interface only and permits arbitrary constitutive laws for the underlying solid, see examples in Shvarts (2019).

2.2 Thin fluid flow

The thin fluid flow in governed by:

in (4a)
on (4b)
on (4c)

where (4a) is the Reynolds equation for isoviscous incompressible Newtonian fluid (Hamrock et al., 2004), note that the tangential relative motion of the solid walls is not considered here whereas the normal approaching is assumed to be quasi-static, is the fluid pressure field defined on the lubrication surface , which is a projection of the fluid-structure interface on the plane , and the operator in (4a) is defined as . Finally, (4b) are the Dirichlet boundary conditions with a prescribed fluid pressure and (4c) are the Neumann boundary conditions with a prescribed fluid flux , defined on and , respectively ( is the outer normal to ), while the fluid flux is given by:


where is the dynamic viscosity. Note that for each point the thickness of the film is computed as the normal gap of the corresponding point .

2.3 Fluid-structure interface

The equilibrium of the solid and fluid tractions on the fluid-structure interface needs to fulfill the following equation:


where the first right-hand side term is the normal traction due to the hydrostatic pressure, while the second one is the tangential traction due to viscous shear stresses in the fluid that act on the solid’s surface (here, it results from the Poiseuille flow), see (Hamrock et al., 2004) for details. Note that the gradient operator here is defined on the lubrication surface , and the second term in (6) is not exactly perpendicular to the outward normal to the surface . Nevertheless, this slight inconsistency is justified by the requirement of small slopes of the surface geometry for validity of the Reynolds equation in certain applications, such as the fluid flow through fractures (Brown et al., 1995), and is often accepted in elasto-hydrodynamic lubrication problems, see, for example, (Stupkiewicz, 2009; Stupkiewicz et al., 2016).

Moreover, in derivation of the Reynolds equation (4a) the thickness of the fluid film is assumed to be much smaller than other length scales, and therefore, the term corresponding to the tangential traction in (6) is often neglected in lubrication problems, see e.g. (Stupkiewicz, 2018). However, in application to sealing problems, studies of the elasto-hydrodynamic lubrication regime show a noticeable effect of the shear tractions on the seal’s leakage, see (Stupkiewicz and Marciniszyn, 2009). Therefore, for the sake of completeness, we will consider both normal and shear components of the fluid-induced traction in the developed framework.

Note also that we do not need to consider no-slip condition (i.e. zero flow velocity at the fluid-structure boundary), as, for example, in (Farhat et al., 1998), since it is already taken into account in the considered form of the Reynolds equation (4a).

2.4 Trapped fluid zones

The hydrostatic pressure , developed in the -th trapped fluid zone, is applied to the surface of the solid body as the normal traction:


However, the pressure is a priori unknown, and the behaviour of the trapped fluid is not governed by the Reynolds equation. Furthermore, the pressure developed in such trapped fluid pocket can be considerably higher than the pressure in the fluid flow, and therefore a model of a incompressible fluid becomes inaccurate here, see discussions in (Shvarts and Yastrebov, 2018b). Therefore, we will consider the model of a compressible fluid with pressure-dependent bulk modulus , where [Pa] and (dimensionless) are model parameters, which is suitable for fluids typically used in lubrication and sealing applications, see also (Kuznetsov, 1985) for details. According to this model, the pressure of the trapped fluid is a non-linear function of the relative change of its volume:


where is the current volume of fluid in the -th zone,

is the volume of the fluid in this zone at the moment when it was formed, and

is the corresponding initial pressure of this trapped fluid. Therefore, since the behaviour of each trapped fluid zone depends on its own set of initial parameters ( and ), it has to be considered separately from others. As was mentioned above, we assume that the fluid occupies the whole free volume between the contacting surfaces. Accordingly, the volume of the fluid in the -th pool is equal to the volume of the gap between the surface and the rigid flat:


where is the outward normal to the surface , and is the normal to the rigid flat, see Fig. 1.

It is important to note, that the presented problem set-up corresponds to the two-way coupling approach, and therefore not only the displacement field depends on the fluid pressure and vice versa, but also the extent of fluid-flow domain and trapped fluid zones may depend on the morphology of the active contact zone, i.e. on the resolution of contact constraints. Additional effort may be necessary for handling edge effects, e.g. enforcing continuity of surface tractions across boundaries between contact and fluid-flow zones , and also between contact and trapped fluid zones . Below we will discuss in detail our recipes of partitioning the interface and handling these problems. At the same time, the one-way coupling approach, which neglects the action of the fluid pressure on surface of the solid, can also be considered in the present problem statement if equations (6) and (7) are omitted.

3 Variational formulation

Before presenting the numerical framework, we discuss briefly the variational statement of the coupled problem formulated in the previous section. We start by outlining contribution of each sub-problem to the balance of virtual work and then provide the variational formulation of the coupled problem for both one- and two-way coupling approaches.

3.1 Weak formulation of the solid mechanics problem with contact constraints

The variational statement of the solid mechanics problem with contact constraints (2d) is well-known, see e.g. (Kikuchi and Oden, 1988), and consists in finding a function :


where denotes the Sobolev space for vector-valued functions with square integrable derivatives, such that:


where is the space for virtual displacements:


Note that in (11) we used the variation of the normal gap function , related to the virtual displacement as:


where is the normal to the rigid flat.

3.2 Weak statement of the fluid flow sub-problem

In order to obtain the weak form of the fluid-flow problem (4c), we follow the standard approach to elliptic (e.g. steady-state heat conduction) equations, see e.g. (Zienkiewicz and Taylor, 1977), which results in the following statement: find a scalar field :


such that


where is the space for scalar virtual functions:


3.3 Weak form of the fluid-structure interface balance

We may compute the work of the fluid-induced tractions on the surface  (6), corresponding to a virtual displacement , as:


and include it to the balance of the virtual work (11) with the negative sign, since the virtual work of surface tractions has a sign opposite to the one of the work of internal forces. It is also important to bear in mind that the traction vector in (6) and (17) is not known a priori and depends on the fluid pressure , its gradient and the displacement field .

3.4 Virtual work of trapped fluid zones

To capture the effect of the trapped fluid pressure on the solid, we recall the thermodynamic definition of the elementary work done on a system, corresponding to an infinitesimal change of its volume. Following this concept, we compute the virtual work of an -th trapped pool on the surface of the solid as:


where the minus sign is used since an increase of the volume of a trapped pool leads to a decrease of its pressure, and consequently, to a release of the energy of the trapped fluid. Since the volume of the fluid inside a trap is a functional of the displacement field as defined by the integral (9), can be treated as its first variation and computed using the directional derivative:


Therefore, the virtual work of the trapped fluid corresponding to a virtual displacement can be expressed as:


which can now be included into the equation for the virtual work (11), taking into account the contribution of each trapped fluid zone separately.

3.5 Variational formulation of the coupled problem

Combining contributions of the sub-problems outlined above, we provide the variational statement of the coupled problem in the spirit of the monolithic approach (Yang and Laursen, 2009; Stupkiewicz, 2018):

Find and such that:




The weak problem statement (21)-(22) is valid for the two-way coupling approach, when both sub-problems have impact on each other. The one-way coupling approximation for the problem under study can also be considered upon following modifications: (i) omit the fluid-induced tractions on the surface of the solid (22b), (ii) neglect the effect of trapped fluid zones (22c), and (iii) assume rigid solid walls while solving fluid-flow equation (21b). Therefore, in case of one-way coupling, instead of (21), we shall have the following equations:


Note that is still required as an input to compute the normal gap, hence we used the subscript “” for , however, for any given displacement field the fluid sub-problem becomes linear under the one-way coupling approach.

4 Computational framework

In this section we discuss the proposed monolithic finite-element framework for the coupled problem. Similarly to the previous section, we discuss handling of each sub-problem separately, and then we present the monolithic resolution algorithm which combines them. We use here the standard mechanical finite element approach for the numerical simulation of the solid deformation, details of which may be found elsewhere. Note that for brevity the same notations are preserved for discretised entities as were introduced in the continuous problem statement. Furthermore, in order to simplify the discussion and concentrate it on the two-way coupling aspects, we use the small deformations and small rotations assumptions, which is (at least partially) justified by the requirement of small slopes of the roughness profile for the validity of the Reynolds equation in certain applications, see e.g. Brown et al. (1995). Nevertheless, the necessary modifications to take into account large deformations and/or large rotations can be added into presented framework. The development of the proposed scheme was undertaken in the finite-element suite Z-set (Besson and Foerch, 1997; Z-set, 2019).

4.1 Mechanical contact

The solution of the considered coupled problem requires partition of the interface into contact, fluid-flow and, possibly multiple, trapped fluid zones, see (1). In order to make the identification of the interface status self-consistent, it appears natural to associate contact elements with faces of the surface , i.e. adopt the “face-to-rigid-surface” discretization approach, rather than the “node-to-rigid-surface” technique, which associates each contact element with a single node of the surface , see (Wriggers, 2006; Konyukhov and Schweizerhof, 2012; Yastrebov, 2013) for more details.

Furthermore, for each point on the surface

we consider the interpolation of the gap and of the normal traction as, respectively:


where is the nodal gap value, is the nodal value of the contact pressure (treatment of which is method-dependent and will be discussed in detail below), is the shape function associated with the node , and is the total number of nodes on surface , see Fig. 2. Note that the same shape functions are used here for interpolation of geometric gap and surface tractions, however, it is not a necessary condition. We used bilinear shape functions associated with quadrilateral faces of the discretised surface, nonetheless, polynomials of higher order can be utilized, see Puso et al. (2008).

In case of the considered interpolation (24) contact constraints (2b) cannot be satisfied point-wise on the surface . To overcome this inconsistency, we follow the mortar approach (Puso, 2004; Puso and Laursen, 2004) and consider the third condition in (2b) in the integral form over the surface :


Substituting (24) into (25) and considering two first conditions of (2b) at every node of the surface , we obtain the following discrete (nodal) form of the contact conditions:


where is termed as the integral (weighted) gap associated with node and is given by:

Figure 2: Sketch of the interface highlighting contact elements: is the potential contact zone, is the active contact zone (active set), is a face associated to one contact element, is the Lagrange multiplier, which represents contact pressure and is attributed to each node of the surface , is the gap function: on and on .

In order to solve the contact problem with constraints (26), we use the augmented Lagrangian method, which combines benefits of the classic Lagrange multipliers method (exact satisfaction of the constraints) and the penalty method (so-called “active set strategy” is not required), see also (Alart and Curnier, 1991; Cavalieri and Cardona, 2013). Thus, we append Lagrange multipliers to each node of the surface , see Fig. 2, and introduce the following augmented Lagrangian functional:


where is the potential energy of deformed solid, while represents the “potential energy” of the contact and is given by:


where is the so-called augmentation parameter. The following notation of the augmented Lagrange multiplier is introduced: , the sign of which defines the contact state of the node: if the node belongs to the active contact zone, while if the node is not in contact. Once the solution is obtained, values of Lagrange multipliers are equal to respective nodal values of the contact pressure . In formulas (28)-(29) we denote by and vectors of nodal displacements and values of Lagrange multipliers, respectively.

The solution of the contact problem is equivalent to the stationary saddle point of the Lagrangian (28), at which its variation vanishes:


where the virtual work of the internal forces is expressed using the directional derivative:


while the second terms in (30) is equivalent to the virtual work of contact forces, cf. (22a), and the third term in (30) ensures that in the active contact zone. In order to derive the contribution of each contact element to the balance of virtual work and determine the element’s status (independently from the neighbouring elements) we compute the integral gap (27) over the face associated with each contact element:


where is the total number of nodes of the face . The weights are calculated as:


where is the Jacobian of the transformation of the physical coordinates on the surface to the face’s coordinates in the parent space :


where is the position of the -th node of the face. Using the Gauss quadrature rules, the integral in (33) is computed as:


where is the number of Gauss points associated with the face , is the weight coefficient of the -th Gauss point, and are its coordinates in the parent space.

In order to find the contribution of each contact element to the balance of virtual works, we calculate the variation of (29):


where is the displacement vector of the node . Note that in accordance with the infinitesimal strain formulation the Jacobian is not variated. We shall term hereinafter an element as active if at least at one of its nodes , and inactive otherwise.

In order to perform linearisation of the problem, we calculate the second variation of the virtual work :


Finally, the virtual work (36) and its variation (37) could be expressed in a compact form, introducing the residual vector and the tangent matrix of a contact element:


where for brevity we use the notations , and , . Ready-to-implement expressions of the outlined components of and are given in the Appendix A.1.

The residual vector and tangent matrix are updated on each iteration of the Newton-Raphson method and added to the corresponding entries of the global residual vector and tangent matrix. Note that in the frictionless case considered here the tangent matrix of the contact element is symmetric, i.e. .

4.1.1 Post-processing computation of the real contact area

Figure 3: Sketch on the computation of the real contact area. Red circles represent nodes with , black ones correspond to . The dashed area represents the contact area computed by summing up areas of all active elements, see (39). The shaded area is obtained by a refined approach of summing up areas corresponding to nodes with , computed using the nearest Gauss point to node , see (40).

The presented above contact element formulation is sufficient for resolution of the contact constraints (26). However, during the post-processing stage, different methods may be applied to compute the real contact area. A possible straightforward approach is to sum up areas of faces associated with active elements, i.e. the ones that have at least one node with (dashed areas in Fig. 3):


where is the number of nodes of the contact element, is the number of Gauss points associated with the face , is the weight coefficient of the k-th Gauss point, and are its coordinates in the parent space.

However, our study showed that this method of computation of the contact area leads to its significant overestimation. We propose here a more precise approach to computing of the contact area: considering separately each contact element, only if at a node , we add up to the contact area contribution from the Gauss point closest to this node (shaded area in Fig. 3):


where and are the weight coefficient and the position of a Gauss point closest to the node . Note that we assumed here that (the number of element’s nodes where Lagrange multipliers are considered) equals to (the number of Gauss points of the corresponding face). However, it might not be the case if, for example, shape functions of different order are used for interpolation of the geometry and of the contact pressure, cf. (24), see also (Puso et al., 2008), and a different refined approach of the real contact area computation will be required. The comparison of two discussed approaches to real contact area computation will be presented below.

4.2 Thin fluid flow

Next, we discuss the implementation of the weak form of the fluid flow problem, see (21b), into the finite-element framework. In order to concentrate the reader’s attention on the aspects of the coupling, we assume that the prescribed fluid flux on , however, a non-trivial Neumann boundary condition could be included in the present formulation in the standard manner.

Figure 4: Sketch of the interface highlighting fluid flow elements: is the lubrication surface (plane ), where the Reynolds equation is defined, is a face on this surface associated to one fluid flow element, is the fluid pressure DOF, added to each node of the surface of the solid (note that the constant fluid pressure across the film thickness is assumed).

Following the assumption of the constant fluid pressure across the film thickness, we attribute a fluid pressure DOF to each node of the surface and define a finite element for the fluid transport problem for each face of the surface , formed by faces of projected on the rigid flat, see Fig. 4. We use the same interpolation for the gap as in the contact problem, see (24), while for the fluid pressure and the “test” function we also have:


where is the number of nodes which belong to an element. Substituting these expressions into (22d) we obtain the following contribution of one fluid-flow element to the balance of the virtual work:


where are coordinates in the parent space, is the Jacobian matrix defined as:


and is its determinant. The second variation, required for the linearisation of the coupled problem, reads:


Note that the variation of the Jacobian matrix is not considered due to assumptions of small deformations and small rotations.

Introducing the residual vector and the tangent matrix of a fluid flow element (explicitly given in Appendix A.2), we may write:


with and similarly . Note the presented formulation of the fluid-flow element was derived for the two-way coupled problem, however, it is also suitable for the one-way coupling. The only required modification is the assumption of the rigid walls of the solid, according to which the variation of the virtual work with respect to the displacement in (4.2) vanishes and consequently .

4.3 Fluid-structure interface

To include the virtual work of fluid tractions on the solid’s surface (22b) into the numerical framework, we associate a fluid-structure interface element with each face of the surface , see Fig. 5. We use the same interpolation of the fluid pressure and the gap, as in (41) and (24), respectively.

Figure 5: Sketch of the interface highlighting FSI elements: is fluid-structure interface on the surface of the solid, is a face associated to one FSI element, is the fluid pressure at each node of this surface and is the outer normal.

Therefore, the contribution of each fluid-structure interface element to the balance of virtual work reads:


where the gradients of shape functions are computed on the projection of the face on the rigid flat, i.e. in the same sense as in Subsection 4.2. The corresponding Jacobian matrix was defined in (43), and the normal is given by:


Note that here is not , but is computed as in (34). The second variation then takes the following form:


where variations of the Jacobian , the matrix and the normal vector are not considered due to assumptions of small deformations and small rotations.

Finally, the associated virtual work and its variation could be expressed in a compact form using the residual vector and the tangent matrix of a FSI element (details can be found in Appendix A.3):


Note that in case of one-way coupling the action of the fluid pressure on the surface of the solid is neglected, so that the virtual work (22b) vanishes and no contribution of the FSI element is included into the global system, i.e.

4.4 Trapped fluid zones

In order to take into account the effect of pressurized volumes of trapped fluid (22c) on the coupled problem, we follow (Shvarts and Yastrebov, 2018b) and use the nonlinear penalty method to simulate the behaviour of the compressible fluid with a pressure-dependent bulk modulus. We discuss here two possible approaches for implementing this model into the finite-element framework. First, we present a “super-element” formulation for a trapped-fluid element containing all faces of one trapped fluid zone . A possible standard finite-element formulation which computes contributions from each face of the trapped fluid zone separately is also stated, and benefits and drawbacks of these two approaches are briefly discussed.

4.4.1 “Super-element” formulation

In the finite-element framework the volume of the gap (9) can be calculated by the following formula:


where the summation with respect to index is performed over all faces of the surface , and the summation with respect to is over all nodes of the -th face. Thus, we denote by vector of displacements of all nodes on the surface , which shall serve as the DOF vector for the trapped fluid “super-element”. Next, is the normal gap computed for the -th node of the -th face and is the shape function associated with this node; are coordinates in the parent space, and is the Jacobian defined in (34). Finally, is the normal to the -th face, which can be computed in the same way as in (47). The integral in (50) can be calculated using the Gauss quadrature rule as:


where is the number of Gauss points associated with the -th face of the surface , is the weight coefficient of the -th Gauss point, are its coordinates in the parent space.

Therefore, using the expression for the virtual work (20), we may write the residual vector and the tangent matrix for the trapped fluid “super-element” as:


where is the tensor product. Note that we used the small strain formulation and neglected the variation of normals, and therefore the matrix of second derivatives is zero, which simplifies the formulation of .

4.4.2 Standard finite-element formulation

Alternatively to the “super-element” formulation presented above, the standard finite-element formulation can be used, according to which the residual vector and the tangent matrix are constructed using separate contributions from each face of the trapped fluid zone. However, in application to the considered problem the standard approach is bound to certain limitations, which will be discussed below.

Indeed, the volume of the gap (50) can be computed as sum of volumes corresponding to each single face:


where is the vector of displacements of nodes of the -th face only. However, the tangent matrix includes the tensor product , see (4.4.1), and therefore is not sparse and cannot be constructed using the standard assembly process, combining contribution from each face separately.

Nevertheless, the standard finite-element formulation can still be used to handle the trapped fluid model, if the method of the Lagrange multiplier and the penalty method are used simultaneously. In order to show that, following (Abaqus 2018, 2018), we consider the contribution of the trapped fluid to the combined Lagrangian for the coupled problem as:


where is an additional Lagrange multiplier, is computed in the same sense as in (9), while represents the volume of the trapped fluid as a function of its pressure, i.e. the inverse of the constitutive relation (8):


Substituting by in (55), we may express the variation of the term