## Abstract

The mean field fluid dynamics and the friction occurring in the wet sliding contact between inhomogeneous surfaces, characterized by a deterministically repeated pattern of microdefects, are modelled within the Bruggeman effective medium theory. By comparing with the results of an accurate numerical homogenization of the flow equations, and with asymptotic solutions, we discuss the validity of the mean field model and its limitations in relation to the occurrence of clustering and interference effects. Finally, an analytical upgrade to the Bruggeman approach, allowing for inclusion of the clustering effect, is presented and discussed.

## 1. Introduction

The microscale mechanics of interacting surfaces has nowadays received consolidated attention from the research community owing to its major implications for long-standing problems, such as tailored friction and reversible adhesion [1], as well as for fundamental research in the fields of cell adhesion and motility on substrates [2]. Indeed, the structural and material properties of most man- and Nature-made surfaces, usually arranged as an ordered or random hierarchy of defects with respect to ideal characteristics [3], have been shown to play a remarkable role in the formation of the macroscopic properties of a generic contact, a topic of great importance in both Nature and technology [1]. In the case of random distribution of surface defects, most of the multi-scale aspects of the contact mechanics have been recently unified under homogenized contact models, which allow one to effectively capture the physics occurring in the wet or dry interactions of real solids [4–9]. For the particular case of (randomly) fluid-mediated interactions, in [6–9] the authors describe the generic contact interface as a two-dimensional time-varying porous medium, the statistical characteristics of which are determined (along the entire range of roughness length scales) by the fluid–asperity and asperity–asperity mutual interactions [6,7]. At the macroscale, the composition of these local interactions largely affects the mean flow and, therefore, the resulting frictional properties of the contact itself. However, as usually occurs, the practice of tailoring surface properties involves the fabrication of ordered one (or at most two) scale micro- (nano-) structures. This is shown in a number of examples in the literature, from the wet contact mechanics of micro- and macromachine element couplings [10,11] to biomedical devices [12,13] and bioinspired functionalized surfaces [14]. In all such applications, the roughness is therefore produced with a quantized spectral content. As a consequence of the different nature of the roughness and, therefore, of the contact itself, the methodologies developed for randomly rough surfaces cannot be applied in this context.

While, in the near future, an increasing interest in the development of novel theoretical tools for non-random contact mechanics is expected, the (relatively) recently emerged topic of fluid dynamics between narrowed interfaces has already been the subject of a number of key investigations. Indeed, several theoretical studies have been devoted to shedding light on the role of an ordered array of surface microstructures on the effective (i.e. collected or measurable at the macroscopic scale) fluid flow, wall friction and slippage [15,16] properties of the confinement, with particular attention to open- (e.g. microfluidic devices [15]) and closed-channel (e.g. sliding contacts [17]) geometries.

In figure 1^{1} the open- and closed-channel configurations are schematically described for a texture made of pillars arranged in a square lattice. In the latter case (figure 1*b*), which is of closer interest for the present work, the investigation is usually focused on determining the optimal characteristics of (only, up to now) structural texture that minimize the frictional dissipation of the contact pair. Note that in such applications, apart from a few cases, the fluid dynamics occurring from the micro- to the macroscale has to be entirely solved in order to predict simple scalar quantities such as the frictional resistance. From a practical point of view, a full-scale numerical resolution, even in the case of thin film flow, can require the development of ad hoc numerical strategies in order to reduce the (high) computational cost [20,21]. Nevertheless, the larger drawback of full-scale approaches basically resides, in my opinion, in their intrinsic lack of understanding of the mean flow dynamics occurring at the interface (see, for example, the discussion reported in note B of [22]).

An alternative to full-scale numerics is given by the homogenization of the contact equations, which would result in a small scale-free contact problem. This possibility has been recently discussed in [23], where the author presented, for the first time, a mean field model of ordered texture hydrodynamics based on the application of the Bruggeman effective medium (BEM) [24,25] approach to thin film flow (in the lubrication approximation). This theory, based on the formal equivalence between the reduced fluid flow formulation and the diffusion problem, allows one to analytically evaluate the average hydrodynamics of the fluid confined at the interface of textured solids in steady sliding. In particular, the homogenized effect of the local texture characteristics (figure 2), such as hole depth and slip length, is quantified in terms of correction factors (known as flow factors; see the electronic supplementary material, S1) to be embedded in an effective thin film flow equation. The model has already been tested in successful comparative studies with both numerical [23,26] and experimental [26,27] results from independent investigations. However, there exist no theoretical arguments supporting the applicability of the BEM approach (which is well known to hold for defects that are randomly distributed in the conducting medium) to an ordered distribution of defects, in our case a lattice of flow conductivity defects. Unravelling the physics behind this would be very important not only for the understanding of the sliding friction for textured surfaces, but also for the generic conduction problem of deterministically distributed phases in a conducting medium.

Therefore, in this work, an attempt is made to elucidate the main mechanisms behind the texture hydrodynamics theory, with the purpose of evaluating the applicability of the BEM features to the diffusion reformulation of the lubrication flow. This is provided by developing a Green’s function approach which allows us to accurately (numerically) homogenize the flow equations at a convenient computational cost. While the numerical tool can be applied to any given texture characteristics for the resolution of the homogenization problem, in this work its use will be restricted to the comparison with the predictions of the BEM texture hydrodynamics (BTH) theory [23]. The outline of the paper is as follows. In §3, after a concise summary of the BTH theory, the numerical homogenization model is described (whose details are reported in the electronic supplementary material, S2). In §4 the numerical homogenization model is used in a step-by-step validation of the BTH assumptions, whereas in §5 an analytically corrected BEM texture hydrodynamics (cBTH) theory, which allows one to obtain a numerically exact homogenization of the flow even at large values of texture area coverage, is presented. Finally, in §6 a discussion on the finite size effects and their influence on texture hydrodynamics predictions is provided.

## 2. The lubrication theory of an array of surface defects: mean field theory

Here the case of a lubricated contact in steady sliding, occurring in the isoviscous rigid regime [28], is considered. In figure 3*a*, the schematic description of a possible contact pair (a sector-pad thrust bearing [29]) is shown. In particular, the latter presents three length scales at which the interaction can be described: (i) a macroscopic scale (represented by a lateral size of the pad) where the mean field theory is solved, (ii) the scale of a representative subset of defects (or inclusions, as they are called in the Bruggeman theory), where the lubrication problem is homogenized, and (iii) the length scale of the generic defect. The textured area is locally (i.e. over the size covered by a representative subset of defects) characterized by an ordered and locally homogeneous lattice of defects (see figure 3 for a schema).

For the previous system, within the lubrication assumption and considering a slower variation of the fluid velocity with the coordinates (figure 1) *x* and *y* as compared with the variation in the orthogonal direction *z* (assuming also a slow time dependence of the Navier–Stokes equations, and neglecting the local effect of cavitation^{2}), it is very easy to show that the pointwise flow *J*(**x**), describing the fluid flow at a generic position **x**=(*x*,*y*) in the contact domain (and including Navier’s boundary condition of slip length *l*(**x**)), reads
*u*(**x**) is the separation, *η*(**x**) is the fluid viscosity, *p*(**x**) is the fluid pressure and **v**_{0} is the sliding velocity (see, for example, the nomenclature in table 1). By collecting prefactors, equation (2.1) can be rewritten as
*Λ*(**x**) and *Γ*(**x**) are, respectively, the generic hydraulic and sliding pointwise conductivities of the interface, expressed in a tensorial form in order to preserve the generality of the model. Equation (2.2), together with the mass conservation equation

In a classical full-scale calculation to determine, for example, the sliding friction, equations (2.1)–(2.3) have to be solved numerically from the macroscale of the contact to the microscale of the surface defects. One way to avoid this is to integrate out the contributions from the surface texture (occurring at the scale of the defects), by applying the two-dimensional BEM approach [24,30] to equations (2.2) and (2.3). This results in the BTH theory [23] summarized in the following. The basic idea is illustrated in figure 3*a*. Owing to the large separation in length scales between the macroscale, characterized by the lateral size *L*_{0} of the pad, and the size *s* of the generic defect (with *L*_{0}≫*s*), the lubricated contact can be conveniently described at the macroscale with a mean field (i.e. slow varying) dynamics, whose effective conductivities are *Λ*_{eff}(**x**) and *Γ*_{eff}(**x**) (with |**x**|∼*L*_{0}). These are related to a reduced set of parameters describing the defects’ characteristics and their distribution along the textured surface (i.e. to *Λ*(**x**) and *Γ*(**x**), with |**x**|∼*s*), and have to be determined by solving the corresponding lubrication problem at the REV scale.

The basic physics behind the BEM approach, from which the REV scale problem is determined, is shown in figure 3*b*. A generic subset of inclusions (i.e. repeated defects or inhomogeneities in the lubrication conductivities) of *constant* conductivities (*Λ*_{i},*Γ*_{i}) is immersed (in the average sense) in the effective medium conductivities (*Λ*_{eff},*Γ*_{eff}). Therefore, the average (effective medium) flow in the subset area, caused by a source at infinity (∇*f*_{0}), will be simply given by the (weighted) superposition of the flows locally induced by the embedded variations of conductivities (i.e. the defects), plus the base conductivity flow. The latter is classically associated with the same effective medium flow, hence resulting in a set of mean field equations in the effective conductivities (see later). It is stressed that, in the original BEM approach, the weighted superposition is obtained through the probability density function describing the distribution of defects, whereas for a deterministic pattern of defects it simply coincides with the summation weighted by the defects area coverage. The main task in identifying the REV problem is, therefore, to analytically determine the flow associated with the defect (of generic in-plane elliptical shape), which is assumed to be *insulated* and *embedded in the effective medium conductivities*.

While the algebra describing the previous REV problem is not really relevant here (which can be found in [23,31]), it is observed that the two fundamental BEM features, whose applicability to our texture problem will be discussed in this work, are the dilution and the embedded assumptions. In particular, the dilute limit, i.e. the neglect of high-order (proximity) interactions between defects (the so-called clustering or aggregation [32]), is theoretically (usually) identified by the constant reduced pressure gradient inside the defect, an approximation which is equivalent for the electrostatic problem to the constant polarization inside the inclusion [31], and which allows one to obtain an analytical resolution of the REV problem. Interestingly, it will be shown, in the following, that the classical dilute assumption not only concerns the constant polarization of the inclusion, but also implicitly admits the occurrence of disruptive interfering flow interactions between defects.

However, by admitting the validity of the dilution and the embedded-inclusion assumptions, as per the classical BEM approach, the REV problem allows one to obtain the following tensorial equations [23] in the effective conductivities *Λ*_{eff} and *Γ*_{eff}:
*Λ*=*Λ*_{i}−*Λ*_{eff}, d*Γ*=*Γ*_{i}−*Γ*_{eff} and **E**_{0} is defined by the geometry of the generic inclusion of index-. **I** is the identity matrix. Finally, the effective Reynolds equation [23], to be numerically solved at the scale of the macroscopic contact, reads
*U*_{m}=*v*_{0}/2 and *p*_{hi} is the -indexed inclusion area density, the latter defined as the coverage area of the generic repeated inclusion with respect to a nominal area. The simplest and commonly encountered case is the two-component texture; see, for example, the lattice cell in figure 4*a* and a (possible) corresponding three-dimensional representation in figure 1*b*. In such a case, the REV problem simplifies to a two-component problem, with *i*=1, 2 in equations (2.4) and (2.5). In particular, one component describes the generic repeated circular geometry (associated with the conductivities (*Λ*_{h},*Γ*_{h})), whereas the other represents the repeated unstructured substrate (associated with (*Λ*_{0},*Γ*_{0})); see figure 4*a* for a schema. Furthermore, *Γ*_{h}=**I***β*_{h}*h*_{h}, where *η*_{h} is the local lubricant viscosity, *h*_{h} is the inclusion height, *l*_{h} is the local Navier slip length, **I** is the identity matrix and [23]

## 3. The lubrication theory of an array of surface defects: numerical averaging model

### (a) General theory for arbitrary defect conductivities

To allow comparison with the predictions of equations (2.4) and (2.5), (2.2) and (2.3) are numerically homogenized.^{3} Let us restrict our attention to the case where the surface texture is constituted by a two-component texture (figure 4*b*). The substrate conductivities (*Λ*_{0},*Γ*_{0}) are assumed constant and principally valued in the reference medium, resulting in a pointwise flow given by *J*_{0}(**x**)=−*Λ*_{0}∇*p*(**x**)+*Γ*_{0}*U*_{m} in *Ω*_{0} (untextured domain; see figure 4*b*). The pointwise flow *J*_{e}(**x**) over the array of surface defects (domain *Λ*_{e}(**x**) and *Γ*_{e}(**x**), reads
*δJ*_{e}(**x**) can be defined by
*J*_{0}(**x**)+*δJ*_{e}(**x**)]=0, i.e.
*b*), by adopting a Green’s function formulation as follows. By first solving ∇⋅[*Λ*_{0}∇*p*(**x**)]=*δ*(**x**), one obtains the following kernel (Green’s function) *K*(**x**):
_{01} and λ_{02} are the diagonal elements of the substrate hydraulic conductivity *Λ*_{0} (note that ∇⋅(*Λ*_{0}∇*K*)=*δ*(**x**)). Finally, the pressure gradient field in the REV, characterized by an *R*_{0}-sized domain of ordered inclusions, can be obtained by solving equation (3.1) with the boundary condition ∇*p*_{0} at *R*_{0}), in
*R*_{0} is large enough with respect to the inclusion representative size, each inclusion far away from the REV borders will be characterized by the same pressure gradient field (owing to translational invariance). It is clear that, in the proximity of the REV borders, the defects pressure gradients will be different from the pressure gradients in the inside population of the cluster. However, this finite size effect, affecting the calculation of the average flow in the REV, introduces an error *ε*∼*N*^{−1}, where *N* is the number of defects per REV side (i.e. assuming a constant lattice size, for

Finally, equation (3.3) can be rephrased with respect to a central inclusion (of domain; *Ω*_{e0}; figure 4*b*) of the REV, resulting in the following Fredholm equation to be solved *on the single lattice cell*:
*p*_{e} is the pressure gradient in *Ω*_{e0}, and the effective kernel **x**_{ci} corresponds to the reference centre of the generic inclusion of the array (figure 4*b*). By adopting the variable change ∇*p*_{e}(**x**)=∇*f*_{e}(**x**) *f*_{e}(**x**) (which has to be solved by a simple matrix inversion for calculating the REV-scale problem):
*Ω*_{1−e0}), the pressure gradient can then be easily calculated as a function of the known ∇*f*_{e}(**x**), by substituting the latter in equation (3.4),
*Ω*_{1}=*Ω*_{1−e0}+*Ω*_{e0}, the average flow can then be calculated *Λ*_{eff} and *Γ*_{eff} (note: 〈*J*〉_{|∇p0|=0}=*Γ*_{eff}*U*_{m}, whereas 〈*J*〉_{|Um|=0}=−*Λ*_{eff}〈∇*p*〉, where

### (b) General theory for defects with constant conductivities

In the easiest case where the conductivities of the array of defects are constant, resulting in *f*(**x**), the hydraulic conductivity *Λ*_{eff} can be determined from (considering that ∇*f*(**x**)=∇*p*(**x**)):
*Γ*_{eff} can be calculated (considering ∇*p*(**x**)=∇*f*(**x**)−∇*f*_{0} and *Λ*_{eff}*Ω*_{1}[〈∇*f*〉−∇*f*_{0}]=*Λ*_{eff}*Ω*_{1}〈∇*p*〉. Note that by adopting two independent boundary conditions (i.e. orthogonal), e.g. ∇*f*_{0} and *f*(**x**) are thus independent (owing to the linearity of the integral operator), allowing for the resolution of equations (3.10) and (3.11) in the effective medium conductivities *Λ*_{eff} and *Γ*_{eff}.

The following sections focus on the comparison between the results of (3.10) and (3.11), and the corresponding BEM equations (2.4) and (2.5). Observe that, because the calculation of both effective conductivity tensors depends exclusively on the reduced pressure gradient, comparative results will be produced, for the sake of brevity, only for the effective hydraulic conductivity *Λ*_{eff}.

## 4. On the Bruggeman effective medium approach applied to texture hydrodynamics

In this section, a discussion will be provided on the physics behind the BEM approach, in particular on its applicability to the wet sliding contact of ordered texturing. As anticipated, the three key BEM assumptions are (i) the neglect of high-order interactions between defects (e.g. the aggregation or clustering effects), (ii) the constant polarization inside the inclusions (i.e. the constant reduced pressure gradient inside the defects), and (iii) the null total polarization of the superposed REV conductivities with respect to the effective conductivities (i.e. the embedded assumption). In the literature of inhomogeneous material conductivity, the aggregation effect, resulting from a certain clustering/chaining of the inclusions, can be taken into account in the BEM, in particular by updating the depolarization factors adopted in the mean field conductivity equations [32,33], equations (2.4) and (2.5) in our case. Here, the equivalent of such depolarization factors for our lubrication problem will be derived, and it will be analytically shown, supported by the numerical validation, that the mutual (high-order) interaction between surface defects is not generally negligible, even for large spacing lattices. This interaction, which could be used to modulate flow and friction (as shown in the cBTH theory presented in §5), is surprisingly exactly zero for the most common lattices of defects.

### (a) The aggregation effect of two interacting circular superficial defects embedded in a homogeneous medium

The high-order mutual interactions between defects (aggregation effect) are neglected in the BTH, i.e. equation (3.8) is analytically solved in the dilute limit, where it is imposed,

In order to test the range of accuracy of such an assumption, equation (3.8) is *analytically* solved without the dilute limit, but keeping the constant-reduced pressure gradient inside the inclusion for the coupling represented in figure 5. To simplify the derivation, the case of a surface defect characterized by a circular shape, with base conductivities λ_{01}=λ_{02}=λ_{0}, is considered. In such a case, equation (3.8) can be rephrased in
**x**_{ci}|=*r* and *E*_{xx}=*E*_{yy}=(2λ_{0})^{−1} and *E*_{xy}=0 if *r*=0, where *r* is the radial distance between the inclusion centres (normalized with respect to the circular inclusion radius), and *θ* is the angular shift; see figure 5 for the schematic, and figure 6 for the contour plots of the depolarization factors (for λ_{0}=1). Observe first that, for a pair of inclusions placed at a distance *r*, the aggregation effect, which can be quantified by the depolarization factors, weighs with *r*^{−2}. Therefore, the mutual influence *cannot be neglected a priori, and not only for closely packed inclusions*. However, interestingly, the periodic dependence of the depolarization factors on the angular offset *θ* may determine a null-resultant mutual interference between inclusions when assembled to form a lattice. As an example, for an array of circular defects arranged in a square lattice, the aggregation effect is remarkably null, i.e. *disruptive interference*, remarkably, explains the success of the BTH in accurately predicting the average flow dynamics of common surface textures [23,26], for example, represented in figure 1*b*.

### (b) The polarization of two interacting circular superficial defects embedded in a homogeneous medium

In order to validate the constant polarization assumption, one can numerically solve equation (3.8) for the previous system of two circular inclusions (figure 5), for λ_{0}=1 and for different values of *r* within the range (with a logarithmic step) [2.01,2.75,…,13.0,17.7], ideally corresponding to a range of texture area densities [0.775,0.416,…,0.0186,0.010]. In particular, in figures 7 and 8 the average (numerically calculated) reduced pressure gradient 〈∇*f*_{e}〉={〈∇_{x}*f*_{e}〉,〈∇_{y}*f*_{e}〉} is shown as a function of the angular offset *θ*, compared with those analytically determined by adopting equation (4.1) for different values of *f*_{e} reads
*is fully justified* in the description of the average reduced pressure gradient inside the inclusions, even at large values of texture area density. As an example, at *p*_{h}=0.416, corresponding to *r*=2.75, there is almost no difference between the analytical and the numerical results. Only for closely packed inclusions (corresponding to the smallest distance in figures 7 and 8), the constant polarization does not strictly hold, as was expected. However, the comparisons are still quantitatively accurate for packed clusters, with an accuracy that increases for defects of smaller magnitudes (

### (c) The mean field flow dynamics of an infinite array of circular superficial defects

In this section, equation (3.8) is numerically solved for a square lattice of circular inclusions (figure 4*b*) with the purpose of comparing the numerically calculated effective hydraulic conductivity with the predictions of equation (2.4) (BTH). Because the square lattice geometry eliminates the aggregation effect (as shown before), this comparison will allow us to check for the last BEM assumption adopted in the BTH, i.e. the ideal representation of the generic inclusion as immersed in the effective medium conductivity (the embedded assumption), which is expected to accurately hold only for randomly and homogeneously distributed inclusions.

For the texture considered in figure 4*b*, the effective hydraulic conductivity is isotropic, with the main diagonal element denoted by λ_{eff}. The latter is obtained from equation (3.10) after numerically solving the REV problem given by equation (3.8). Therefore, in figure 9a,*b* the dimensionless hydraulic conductivity (known as the pressure flow factor) λ_{eff}/λ_{0} is shown as a function of the texture area density *p*_{h} (dots, from equation (3.8)), compared with the BTH (solid line) and asymptotic (dashed line) predictions (the asymptotic model is derived in the electronic supplementary material, S3). Different figures are for different values of _{0}=1), in particular *a*) and *b*) (note that

It is interesting to observe that, although the analytically and numerically calculated pressure gradients (figures 7 and 8) match to a high accuracy, the BTH flow factors (figure 9*a*,*b*) do not show the same remarkable agreement. While this could be an inaccuracy expressed by the numerical tool at increasing *p*_{h}, this is not the case here, because in figure 9*a* (*b* where, however, differences from the numerical and the analytical solution are smaller in magnitude.

Therefore (despite both cases showing physically acceptable differences between the numerical and analytical flow factors), the error introduced in the flow factors is not related to the accuracy of the analytical expressions of the average pressure gradient adopted in the BTH, but is undoubtedly a consequence of the third BEM assumption (the embedded assumption; see above). Indeed, the original BEM approach was intended for homogeneously *and randomly* distributed inclusions where, therefore, a generic inclusion is expected to be (averagely) immersed in the (average) effective medium conductivity. In the case where defects are only homogeneously distributed, as in the case of texture hydrodynamics, the defects, even in the case of multi-modal arrangements, are always surrounded by the base conductivities (unless in the case of closely packed defects). As might be expected, this should affect the average flow, with respect to the BEM assumption, especially for increasing values of texture area density. Therefore, it could be very useful to correct the BTH in order to take into account such a phenomenon. As will be shown, the cBTH theory can largely improve (in numerical accuracy) the flow factors’ calculations with respect to the original theory; however, its use is practically justified only for relatively large values of defect magnitudes (i.e. for large |dλ|), where the BTH starts to deviate from the numerical predictions.

## 5. A corrected formulation of the texture hydrodynamics theory (cBTH)

Here, the notations *Λ*_{0},*Γ*_{0}) and (*Λ*_{i},*Γ*_{i}) are, respectively, the base and the generic defect tensorial conductivities. Thus, the effective flow *J*_{eff} acting in the REV can be determined from
*In the case of a two-component texture* (see, for example, figure 4*a*), by substituting equation (5.2) into equation (5.1), one obtains two tensorial equations for the two unknown effective conductivities *Λ*_{eff} and *Γ*_{eff}:
**E**_{0eff} will be isotropic too, i.e. **E**_{0eff}=(2λ_{0})^{−1}**I**, and *Λ*_{eff}=λ_{eff}**I**, *Λ*_{e}=λ_{e}**I**. Substituting the latter in equation (5.3), after some manipulation one obtains
_{eff}/λ_{0} is shown as a function of the texture area density *p*_{h} for *a*) and *b*). Note that the cBTH predicts results in very good quantitative agreement with the numerical ones, even to large values of the texture area density.

Furthermore, by using equations (5.3) and (5.4), one can accurately calculate the effective conductivities in the case where the aggregation effect is non-negligible, such as for circular defects embedded to form a rectangular lattice of aspect ratio *γ*=*l*_{x}/*l*_{y}; see figure 11 for the schematic. In such a case, *E*_{0eff} has to be calculated with a Peklenik number *γ* [23], for which one gets *E*_{0eff,xy}=0 and:
*e*_{xx,γ}=(1+*γ*)^{−1} and *e*_{xx,γ}=(1+*γ*^{−1})^{−1}. Finally, by using equation (5.3), one obtains two equations for the effective conductivities λ_{eff,x} and λ_{eff,y}:

In figures 12 and 13 the main diagonal elements of the pressure flow tensor λ_{eff,x}/λ_{0} and λ_{eff,y}/λ_{0} (figures 12*a* and 13*a*), as well as the corresponding average reduced pressure gradients within the inclusions 〈∇*f*_{e}〉 (figures 12*b* and 13*b*), are shown as a function of the lattice aspect ratio *l*_{y}/*l*_{x} (see figure 11). Results are from cBTH (equation (5.6)) and equation (3.10), with λ_{0}=1. Moreover, in each displayed set of calculations four values of texture area density are shown, in particular *p*_{h}=0.375, 0.192, 0.098, 0.050. (Note that larger values of *p*_{h} correspond to a smaller range of variation of the lattice aspect ratio *l*_{y}/*l*_{x}. This is not indicated in the figures for the sake of brevity.) Observe that, apart from very small and very large *l*_{y}/*l*_{x} values, where the inclusions are closely packed, the cBTH analytical results are in excellent agreement with the numerical ones.

Finally, the locally averaged wall shear stress, acting on the untextured surface, reads
*Φ*_{fp,eff} and *Φ*_{fs,eff} can be calculated.

## 6. Discussion

The BTH theory has been successfully tested in a number of comparisons with both numerical [23,26,27] and experimental [26] results, showing, in particular, that both friction as well as supported load can be accurately calculated within the model. For a single thrust bearing pad, partially textured with square holes arranged in a square lattice, in figure 14 the cross-section pressure field is shown as a function of the contact position (along the sliding axis) as obtained from independent exact numerical calculations [20] (red dotted curves), and is compared with the predictions (with same texture area density but for circular holes) for the BTH [23] (dashed black line) and cBTH (solid line) theories. Note that the cBTH curve is characterized by smaller pressure values with respect to the BTH predictions (as must be expected; see the following); however, these are closer to the numerical results.

In the case of striped defects (e.g. structural grooves, slippery and structural groove, etc.), the cBTH formulation is equivalent to the original theory (the derivation is omitted here for simplicity). Moreover, in such a limiting case, no aggregation effect occurs, so that both theories make use of the correct depolarization factors for the dilute limit. In the remaining part of this section, the BTH for striped defects will be used to shed light on a particular and poorly understood aspect of texture hydrodynamics, i.e. the role of the representative in-plane size *s* of the generic defect. As shown in this work, as well as in [23,26,27], the in-plane size of the generic defect is not relevant to the flow factor calculations, because in the thermodynamic limit there cannot be any dependence scaling with *s*/*R*_{0} (i.e. a finite size dependence), where *R*_{0} is the size of the REV. However, clearly, an upper length scale describing the contact interface is always questioning the range of application of any mean field theory. This sometimes leads to inaccurate considerations about the validity of the same mean field theories, as a consequence of inappropriately inferring general theoretical conclusions from numerical calculations (see, for example, the discussion in [34] for the case of random roughness contact mechanics).

Here, the finite size effect is quantified in the case of striped defects due to the availability, at least for the lubrication of transverse grooves, of exact analytical results [35]. In particular, in [35] the authors report on the lubrication model for a slider bearing (with parallel nominal surfaces) which is partially textured with transverse grooves. The one-dimensionality of the flow and the adoption of a stepped texture profile allows for the piecewise analytical integration of the Reynolds equation and, therefore, for the determination of the macroscopic quantities (as friction or supported load) in terms of summation over the number of profile steps. Thus, in order to have an estimation of the finite size effect, a comparison with the predictions of such a theory is particularly useful. In figure 15*a* the dimensionless load support *p*_{h} and groove ratio *h*_{h}/*h*_{f}. *F*/*B* is the load by transversal length, *h*_{min} is the (nominal or) minimum fluid gap thickness, *η* is the dynamic viscosity, *v*_{0} is the sliding velocity and *L* is the bearing length in the sliding direction. In the figure, the horizontal (asymptotic) lines are given by the BTH theory, whereas the remaining are calculated in [35]. Very interestingly, at increasing number of grooves, corresponding to a decrease in the groove in-plane size *s*, all load curves converge to the BTH limiting values. Inversely, for larger groove sizes, the analytical result [35] diverges from the BTH and correctly predicts the load support of the Rayleigh step, which is well known to be the optimal structural solution that maximizes the supported load for such macrogeometry of the contact. This result suggests two major outcomes. Indeed, if, on the one hand, the finite size effect increases the supported load, i.e. optimizing with the BTH is a safe procedure, then, on the other hand, increasing the plain bearing load support by means of structural microtexturing is not as effective as Rayleigh’s (two-dimensional) or Akers’ (three-dimensional) macrogeometry. While, apparently, the latter considerations may question the proficiency of texturing in practical applications, it is stressed that the microstructural modification can have a key role in increasing the contact performances (i) when intended as a targeted (locally optimized) modification of a given nominal contact geometry or (ii) when combined with other sources of flow conditioning, as a slippage microtexture [26,27]). In the latter case, which is widely discussed in [27], the load support can be much larger than any micro- or macrostructural modification could allow.

It is also interesting to determine the friction force corresponding to the contact geometries reported in figure 15*a* and, in particular, the dimensionless friction force *α* is the texture area coverage), i.e. the friction force is not sensitive to the in-plane size of the groove and, therefore, no finite size effect occurs. In terms of the friction coefficient, in figure 15*b,* *p*_{h} and for a fixed groove ratio and number of grooves. As a consequence of the previous results, the friction coefficient, as well as its minimum as a function of the contact characteristics, will be accurately determined by the BTH theory, which therefore is confirmed to be an effective tool for the fast and accurate optimization of frictional surface textures.

## 7. Conclusion

I have discussed the physics behind a recently developed mean field theory of texture hydrodynamics, based on the BEM approach. The core BEM assumptions, of relevance for the lubrication problem, have been successfully tested, and the range of applicability of the model has been definitely assessed, both with ad hoc numerical calculations and by comparing with the results of independent analytical investigations. Moreover, I have also introduced a novel numerical approach for the fast and accurate homogenization of texture-induced flow based on Green’s function formulation of the lubrication problem. Finally, analytical corrections for the original texture hydrodynamics theory, which allow one to take into account aggregation effects as well as to obtain numerically accurate flow factors even to large values of texture area density, have been presented and discussed.

## Data accessibility

Supplementary information is available online. Correspondence and requests for materials should be addressed to M. Scaraggi (michele.scaraggi{at}unisalento.it).

## Competing interests

I have no competing interests.

## Funding

The research has not been supported by external funding.

## Acknowledgements

The author sincerely acknowledges the useful remarks and the manuscript text revision as kindly performed by the referees. The author also acknowledges his wife’s immense patience for having partially unloaded him from newborn duties, resulting in some free time for this publication.

## Footnotes

↵1 In this work, the structural defect is assumed rigid (see also the following), i.e. the defect deformation induced by the same fluid–structure interaction is neglected. This may be the case of (locally) isoviscous rigid interactions (as in plain bearings); however, for soft contacts (as typically occurs in biological interfaces, e.g. in the lubrication of polymer brushes [18]) the elastic coupling cannot usually be neglected. A stimulating discussion on the inclusion of pillar–elasticity effects on the calculation of hydraulic conductivities (poroelasticity) for wet carpets and brushes has been recently presented in [18]. However, the latter theory [18] cannot be applied here owing to the different microgeometry (and microstructure aspect ratio) involved in the present study which, moreover, requires the effective Couette term to be calculated. Furthermore, microelastohydrodynamics calculations (not shown here for brevity) with pillar microstructures show that soft pillars may deform superficially assuming a wedge configuration, similar to a tapered microbearing pad. While this is not expected to affect the effective flow conductivities, the local shear stresses may experience a substantial variation, as typically encountered in microelastohydrodynamics [19]. Further investigation is clearly required on this topic.

↵2 The generic inclusion height which can be accurately calculated within the model is limited to the applicability of the thin flow assumption (lubrication approximation). Therefore, microstructured surfaces produced, for example, with pulsed laser ablation fit well with the theory, whereas polymer brushes, usually characterized by high aspect ratio values, cannot be described within the thin flow assumptions. A quantitative discussion on the effect of the aspect ratio on the flow and stress factors can be found in [26], §7.2. Moreover, in the same paper [26], §7.1, it has been demonstrated that cavitation, if occurring only at the scale of the defect, marginally affects the flow factor curves, but it substantially alters the shear stresses (owing to the proportionality of the average viscosity to the average density of the cavitation mixture). Therefore, cavitation when occurring locally can be accounted for. The case when cavitation propagates to the macroscale (through, for example, the initial formation of cavitation fingers) is different. However, for such a case, no study is yet available demonstrating the suitability of the BTH theory. Further investigation on this topic is needed.

↵3 It is stressed that the numerical procedure to be discussed is of general applicability, allowing one to solve the REV-scale problem for any texture characteristics, such as non-constant defect conductivities (as occurring, for example, in substructured dimples, or simply where the dimple is a paraboloid, etc.), multi-modal array of defects (i.e. textures constituted by an ordered repetition of differently shaped defects), etc. This extended topic, however, is beyond the purpose of the present study, and it will be the subject of a dedicated contribution.

- Received September 30, 2014.
- Accepted October 12, 2015.

- © 2015 The Author(s)