Analysis of the thin film flow in a rough domain filled with micropolar fluid

Computers and Mathematics with Applications 68 (2014) 1915–1932 Contents lists available at ScienceDirect Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa Analysis of the thin film flow in a rough domain filled with micropolar fluid Igor Pažanina,∗, Francisco Javier Suárez-Graub aDepartment of Mathematics, Faculty of Science, University of Zagreb, Bijenička 30, 10000 Zagreb, Croatia bDepartamento de Ecuaciones Diferenciales y Análisis Numérico, Universidad de Sevilla, C/ Tarfia s/n, 41012 Sevilla, Spain article info Article history: Received 24 March 2014 Received in revised form 2 October 2014 Accepted 3 October 2014 Available online 18 October 2014 Keywords: Thin-film flow Micropolar fluid Rough boundary Different scales Asymptotic expansion Two-scale convergence abstract Inspired by the lubrication framework, in this paper a micropolar fluid flow through a rough thin domain is studied. The domain’s thickness is considered as the small parameter ε, while the roughness is defined by a periodical function with period of order ε2. Starting from three-dimensional micropolar equations and using asymptotic analysis with respect to ε, we formally derive the macroscopic model clearly detecting the effects of the specific rugosity profile and fluid microstructure. We provide the rigorous justification of our formally obtained asymptotic model by deriving the effective system by means of the two-scale convergence. ©2014 Elsevier Ltd. All rights reserved. 1. Introduction The classical lubrication problem is mainly concerned with the situation in which two solid surfaces being in relative motion are separated by a thin layer of fluid acting as a lubricant. Such situation appears naturally in applications consisting of moving machine parts, namely the journal bearings. Fluid film bearings are machine elements whose function is to promote smooth relative motion between two surfaces and are crucial factors in limiting the dissipation of energy. The ultimate goal is that fluid film bearing is well designed so that the wear is not an issue (two surfaces are completely separated by the lubricant). For that reason, it is essential to understand the behavior of the fluid film in such machine elements. The first result goes back to Reynolds and his celebrated work [1] published in 1886. He studied the thin film flow in a rather heuristic manner and did not provide any relation between his model and the Navier–Stokes equations. The formal relationship between Navier–Stokes equations and Reynolds equation in a thin domain was established more than 60 years later in [2,3], while the rigorous mathematical justification of the Reynolds equation for a Newtonian flow between two plain surfaces can be found in [4]. If the gap between the moving surfaces becomes very small, the experimental results from the tribology literature (see e.g. [5–7]) suggest that the fluid’s internal structure should be taken into account as well. A possible way to acknowledge such experimental findings is to employ the micropolar fluid model. Being originally proposed by Eringen [8] in the 60s, the theory of micropolar fluids has gained much attention since it successfully describes the effects of local structure and micro-motions of the fluid elements that cannot be captured by the classical Navier–Stokes model. Physically, micropolar fluids represent fluids consisting of rigid, spherical particles suspended in a viscous medium, where the deformation of fluid ∗Corresponding author. E-mail addresses: [email protected],[email protected] (I. Pažanin), [email protected] (F.J. Suárez-Grau). http://dx.doi.org/10.1016/j.camwa.2014.10.003 0898-1221/©2014 Elsevier Ltd. All rights reserved. 1916 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 particles is ignored. They are, in fact, non-Newtonian fluids with nonsymmetric stress tensor. In view of that, the related mathematical model introduces a new vector field, the angular velocity field of rotation of particles (microrotation) and one new (vector) equation coming from the conservation of the angular momentum. As a result, a complex coupled system of PDEs is obtained, representing a significant generalization of the Navier–Stokes equations. We refer the reader to the monograph [9] (and the references therein) providing a detailed derivation of the micropolar equations from the general constitutive laws together with an extensive review of the mathematical theory and the applications of this particular model. Engineering practice also indicates that it is of great interest to combine the lubrication phenomena with the analysis of the roughness effects. Usually it means that the lower surface is assumed to be perfectly smooth, but the upper is rough and described by a given function. Expressing the boundary roughness using a periodic function, thin-film flow of Newtonian fluid has been extensively studied for different rugosity profiles. The classical assumption is that the size of the roughness is of the same order as the film thickness, i.e. hε(x)=εhx,x ε,0< ε ≪1.(1) In such setting, the effective model turns out to be the classical Reynolds equation (see e.g. [10,11]) and one needs to compute the correctors in order to detect the roughness-induced effects. Same result is obtained for hε(x)=εhx,x εβwith β < 1 (see [12]). In view of that, Bresch and co-authors [13] in 2010 considered a new framework, namely hε(x)=εhx,x ε2.(2) As a result, they derived the asymptotic model in which an extra term (appearing due to the boundary roughness) modifies the standard Reynolds equation at the main order. Whole asymptotic expansion (at any order) of the solution has been rigorously derived in [14] providing the optimality with respect to the truncation error. It is important to emphasize that, roughness pattern described by hε(x)=εhx,x εβwith β > 1 is physically relevant and realistic (see e.g. [15]), and, therefore, has been studied for different situations in recent years. Focusing on the wall laws, the effects of the above setting on the asymptotic behavior of the Navier–Stokes system have been investigated in [16]. Using the asymptotic approximation from [13] derived for the hydrodynamic part of the system, the roughness effects on the heat conduction in a thin film flow have been studied in [17]. A semilinear parabolic problem in a thin rough domain assuming different order to the period of oscillations on the top and the bottom of the boundary has been addressed in [18]. Our goal is to extend the analysis presented in [13] to a case of lubrication with incompressible micropolar fluid. There are not many papers in the existing literature dealing with the mathematical modeling of micropolar fluid film lubrication. Interesting result can be found in [19] where the authors consider a specific slider-type bearing. After writing the governing problem in non-dimensional form, they formally obtain a generalized version of the Reynolds equation in a critical case when one of the non-Newtonian characteristic parameters has specific (small) order of magnitude. Rigorous derivation of such result was brought 14 years later in [20] for two-dimensional setting (see also [21] for micropolar flow in a curved channel). The 3D lubrication problem was recently addressed in [22] and new, second-order Brinkman-type asymptotic model has been proposed. In the above papers, the roughness effects were not taken into account, i.e. the height of the channel is assumed to be of the form hε(x)=εh(x). To our knowledge, the first (and only) rigorous result on the micropolar fluid film lubrication in a thin domain with rough boundary can be found in the recent paper by Boukrouche and Paoli [23]. They consider a micropolar flow in a two-dimensional domain assuming that the height of the channel is given by (1). Employing two-scale convergence technique, they derive the limit problem describing the macroscopic flow. In the present paper, we are going to study a micropolar fluid flow in a three-dimensional domain given by Ωε=(x,z)∈R2×R:x∈ω, 0<z<hε(x),(3) where the height hεis defined by (2). From the point of view of asymptotic analysis, we find this framework more challenging than the classical one (given by (1)) due to the technical difficulties caused by the specific height profile. The main problem related to a fluid flow through a domain with roughness is to deduce in which way the irregular boundaries affect the flow. This is especially important with regard to numerical computations: indeed, roughness is in general too small to be captured by the discretization grid of the simulations. To overcome this difficulty, one can employ the homogenization theory. In view of that, the idea is to replace the irregular domain by a smooth one, and then describe the averaged effect of the roughness in the limit (homogenized) model. For that reason, homogenized models have been of practical interest in numerical codes. In our particular case, starting with original problem (6)–(11) posed in thin rough domain Ωε, we apply the suitable change of variables, namely Z=z/hε(x), to transform Ωεinto Ωwhich is smooth. Then by means of a two-scale convergence technique, we obtain the simplified limit problem posed in Ω, in which the effects of roughness can be clearly observed (see Section 2.3). The paper is organized as follows. After formulating the problem in Section 2, in Section 3we perform a formal asymptotic analysis with respect to the small parameter ε. Introducing a suitable change of variables which takes into account the rough oscillations, we rewrite the governing problem in the ε-independent domain and employ two-scale expansion technique. Since the problem is coupled, we construct the asymptotic expansion of the solution by simultaneously treating boundaryvalue problems for velocity and for microrotation. As a result, we obtain an effective system describing the macroscopic flow and observing clearly the effects of the rugosity profile and fluids microstructure. I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 1917 Fig. 1. The different scales related to the domain. Finally, Section 4is devoted to a rigorous justification of the formally obtained asymptotic model. We apply a convenient variant of the two-scale convergence and verify the effective equations obtained in a formal way. To conclude, we believe that the presented result could be instrumental for understanding the effects of the rough boundary and fluid microstructure on the lubrication process. In view of that, more efficient numerical algorithms could be developed improving, hopefully, the known engineering practice. 2. Formulation of the problem and the statement of the main result 2.1. The domain We consider the fluid flow in the following three-dimensional domain Ωε=(x,z)∈R2×R:x∈ω, 0<z<hε(x).(4) Here we assume that ωis a smooth bounded subset of R2and hε(x)=εh1(x)+ε2h2x ε2.(5) We also define Ω=ω× ⟨0,1⟩ ⊂ R2×R, and denote by T2the torus of dimension 2. As we can see, lower surface is supposed to be plane, while the roughness of the upper surface is described by the given function hε. The functions h1,h2appearing in (5) are assumed to be regular: the positive function h1∈H2(ω) represents the main order part of the roughness, while the T2-periodic function h2∈H2(T2)(with 0 as average in T2) describes the oscillating part (see Fig. 1). 2.2. The equations and boundary conditions In view of the application we want to model, we can assume a small Reynolds number and neglect the inertial terms in the governing equations. Thus, we assume that the flow in Ωεis governed by the following linearized equations: −(ν +νr)1uε+ ∇pε=2νrrot wε,(6) div uε=0,(7) −(ca+cd)1wε−(c0+cd−ca)∇div wε+4νrwε=2νrrot uε.(8) The unknown functions are uε,wεand pεrepresenting the velocity, the microrotation and the pressure of the fluid respectively. Positive constants ν, νr,c0,ca,cdare the viscosity coefficients: νis the usual kinematic Newtonian viscosity, while νr,c0,ca,cdare new viscosities connected with the asymmetry of the stress tensor and, consequently, with the appearance of the microrotation field wε. For the sake of notational simplicity, external forces and moments are neglected and fluid density is assumed to be one. The aim is to study the lubrication process where two rigid surfaces are in relative motion and are separated by a thin layer of fluid. Therefore, we impose the following boundary conditions for the velocity: uε=0 for z=hε,uε=gfor z=0.(9) Here g∈R3is a given constant corresponding to the imposed horizontal velocity of the plane wall. Obviously, g·k=0 implying u3 ε|z=0=0. Here and in the sequel (i,j,k)denotes the standard Cartesian basis. Along the lateral boundary, several types of boundary conditions can be considered, depending on the particular device to be considered. One can use standard Dirichlet boundary condition for the velocity (see [11]), mixed (Dirichlet–Neumann) type condition for the velocity (see [13]) or even combination with pressure boundary condition (see [22]), namely uε×n=0,pε=qεfor x∈∂ω, (10) for given outer pressure qε=ε−2qand normal unit vector n. 1918 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 Finally, to close up the governing problem, we need to prescribe the boundary conditions for the microrotation. Though, recently, some other types of boundary conditions for the microrotation can be found in the mathematical literature (see e.g. [24]), using simple zero boundary condition still seems to be a common practice. Therefore, we impose wε=0 on ∂Ωε,(11) meaning that the fluid microelements cannot rotate on the solid surface. The existence and uniqueness of the solution (uε,pε,wε)∈H1(Ωε)3×L2(Ωε)/R×H1 0(Ωε)3to the described boundaryvalue problem can be established using standard techniques (see e.g. [9,25]). Our goal here is to find the macroscopic law describing the effective flow in Ωεvia asymptotic analysis with respect to the small parameter ε. 2.3. The main result We find the flow to be governed by the following equations    −(ν +νr)∂2 Zv+h2 1∇xp+(ν +νr)MZ∂Zv=0, ∂Zp=0, −(ca+cd)∂2 Zw−2νrh1(−∂Zv2i+∂Zv1j)+(c0+2cd)MZ∂Zw=0 with divxh3 1 12(ν +νr)∇xp=divxh1 2CMg. Here coefficient Mis given by M=T2 |∇Xh2|2dX, while CM=B 6Awith A=1 MeM/21 0 e−Mt2/2dξ−1−1 M(eM/2−1)1 0s 0eM(s2−t2)/2dξds 1 0eMs2/2ds , B=1 M(eM/2−1)1 1 0eMs2/2ds. Note that coefficients Mand CMare both provided in the explicit form and that they depend only on the form of the rugosities. The micropolar nature of the fluid appears through the viscosity ν+νr, while gis the constant corresponding to the imposed horizontal velocity of the plane wall. The above equations have been first formally obtained and then rigorously confirmed via two-scale convergence (see Theorem 1, Section 4). It is important to emphasize that each unknown in the limit problem depends on xand Z, but the pressure depends only on x. For that reason, the above problem can be seen as linear second-order ODE with respect to Zfor vand wleading to explicit expressions for vand w(see Sections 3.5 and 3.6). Since both vand wdepend on the pressure p, it is necessary to deduce an equation for pas well. We proceed in a standard manner and as a result obtain the generalized Reynolds equation posed in ωsatisfied by p(see Section 3.5). Finally, solving the Reynolds equation for pwe can deduce both the velocity vand the microrotation w. To point out how the geometry and roughness of the thin domain affect our problem, let us recall the effective equations in a thin domain without roughness (see [22]). Using the superscript·for the solutions of such problem, the equations read    −(ν +νr)∂2 Z v+h2 1∇x p=0, ∂Z p=0, −(ca+cd)∂2 Z w−2νrh1(−∂Zv2i+∂Zv1j)=0 with divxh3 1 12(ν +νr)∇x p=divxh1 2g. Comparing the above two systems, notice that the roughness of the boundary introduces a new term modifying momentum equations for the velocity and microrotation through the coefficient M, and giving a modified Reynolds equation through the coefficient CM. As in [13], it still remains true the relation p=CM p. Therefore, to our opinion, this contribution represents an important generalization of the results provided in [13,22]. I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 1919 3. Formal asymptotic analysis 3.1. Rescaling As a first step, we need to rewrite the starting problem in the fix (ε-independent) domain. To accomplish that, we first introduce a fast variable X=x ε2capturing the oscillating phenomena of the thin domain. In view of that, the height hε becomes h(x,X)=εh1(x)+ε2h2(X). (12) Next we introduce a new vertical variable Z=z h(x,X)and, correspondingly, the new unknown functions: velocity uε(x,z)= ˜ u(x,X,Z), microrotation wε(x,z)=˜ w(x,X,Z)and pressure pε(x,z)=p(x,X,Z). In the sequel, we also adopt the following notation ˜ u=(v, v3)∈R2×R,˜ w=(w, w3)∈R2×R. The boundary conditions satisfied by velocity and microrotation after performing the change of variables read as follows ˜ v=0 for Z=1,˜ v=gfor Z=0,˜ w=0 for Z=0,1, where g∈R3is given in (9). Moreover, motivated by the periodic nature of h2, we assume that ˜ u,˜ wand pare T2-periodic functions in the variable X, i.e. ˜ u(x,X+1,Z)=˜ u(x,X,Z), ˜ w(x,X+1,Z)=˜ w(x,X,Z), p(x,X+1,Z)=p(x,X,Z). (13) Now we have to express each differential operator appearing in Eqs. (6)–(8) acknowledging the above change of variables. The first and second derivatives of the function θε(x,z)=θ(x,X,Z)can be written as ∇xθε= ∇xθ+1 ε2∇Xθ−1 h∇h·Z∂Zθ, ∂zθε=1 h∂Zθ, ∆xθε=∆xθ+2 ε2∇x· ∇Xθ+1 ε4∆Xθ−2 ε2 ∇h h·Z∇X∂Zθ−1h hZ∂Zθ +|∇h|2 h2Z∂Zθ−2∇h h·Z∇x∂Zθ+|∇h|2 h2Z2∂2 Zθ+1 h2∂2 Zθ, ∂2 zθε(x,z)=1 h2∂2 Zθ, where ∇h(x,X)=ε∇xh1(x)+ ∇Xh2(X)and 1h(x,X)=ε∆xh1(x)+1 ε2∆Xh2(X). (14) The change of variables applied to rotation and divergence yields rot fε=1 ε2rotXf3+(∂X1f2−∂X2f1)k+1 h−∂Zf2i+∂Zf1j+rotxf3+(∂x1f2−∂x2f1)k, div fε=divxf+1 ε2divXf−1 h∇h·Z∂Zf+1 h∂Zf3, for a vector function fε(x,z)=f(x,X,Z)and with rotx(f3)=∂x1f3i+∂x2f3j,rotX(f3)=∂X1f3i+∂X2f3j, divx(f)=∂x1f1+∂x2f2,divX(f)=∂X1f1+∂X2f2. Using the above expressions, we deduce ∇div wε= ∇x(divxw)+1 ε2∇X(divxw)−1 h∇h·Z∂Z(divxw) +1 ε2∇x(divXw)+1 ε4∇X(divXw)−1 ε2h∇h·Z∂Z(divXw) −1 h1hZ∂Zw+|∇h|2 h2Z∂Zw−1 h∇h·Z∇x∂Zw−1 ε2h∇h·Z∇X∂Zw+|∇h|2 h2Z2∂2 Zw −1 h2∇h·∂Zw3+1 h∇x∂Zw3+1 ε2h∇X∂Zw3−1 h2∇h·Z∂2 Zw3 +1 h∂Z(divxw)k+1 ε2h∂Z(divXw)k−1 h∇h·∂Zw k −1 h2∇h·Z∂2 Zw k +1 h2∂2 Zw3k. 1920 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 In view of the preceding calculations, the momentum equation (6) has the following form in an ε-independent domain Ω×R2= {(x,X,Z)∈R2×R2×R:x∈ω, 0<Z<1}: (ν +νr)−h2∆xv−2 ε2h2∇x· ∇Xv−1 ε4h2∆Xv+2 ε2h∇h·Z∇X∂Zv+h1hZ∂Zv− |∇h|2Z∂Zv +2h∇h·Z∇x∂Zv− |∇h|2Z2∂2 Zv−∂2 Zv+h2∇xp+1 ε2h2∇Xp−h∇hZ∂Zp =2νr ε2h2rotXw3+2νrh−∂Zw2i+∂Zw1j+2νrh2rotxw3,(15) (ν +νr)−h2∆xv3−2 ε2h2∇x· ∇Xv3−1 ε4h2∆Xv3+2 ε2h∇h·Z∇X∂Zv3 +h1hZ∂Zv3− |∇h|2Z∂Zv3+2h∇h·Z∇x∂Zv3− |∇h|2Z2∂2 Zv3−∂2 Zv3+h∂Zp =2νr ε2h2∂X1w2−∂X2w1+2νrh2∂x1w2−∂x2w1. The divergence equation (7) in the rescaled domain reads: hdivxv+1 ε2hdivXv− ∇h·Z∂Zv+∂Zv3=0.(16) Finally, the angular momentum equation (8) can be rewritten as follows: (ca+cd)−h2∆xw−2 ε2h2∇x· ∇Xw−1 ε4h2∆Xw+2 ε2h∇h·Z∇X∂Zw +h1hZ∂Zw− |∇h|2Z∂Zw+2h∇h·Z∇x∂Zw− |∇h|2Z2∂2 Zw−∂2 Zw +(c0+cd−ca)−h2∇x(divxw)−1 ε2h2∇X(divxw)+h∇h·Z∂Z(divxw) −1 ε2h2∇x(divXw)−1 ε4h2∇X(divXw)+1 ε2h∇h·Z∂Z(divXw) +h1hZ∂Zw− |∇h|2Z∂Zw+h∇h·Z∇x∂Zw+1 ε2h∇h·Z∇X∂Zw − |∇h|2Z2∂2 Zw+ ∇h·∂Zw3−h∇x∂Zw3−1 ε2h∇X∂Zw3+ ∇h·Z∂2 Zw3+4νrh2w =2νr ε2h2rotXv3+2νrh−∂Zv2i+∂Zv1j+2νrh2rotxv3,(17) (ca+cd)−h2∆xw3−2 ε2h2∇x· ∇Xw3−1 ε4h2∆Xw3+2 ε2h∇h·Z∇X∂Zw3 +h1hZ∂Zw3− |∇h|2Z∂Zw3+2h∇h·Z∇x∂Zw3− |∇h|2Z2∂2 Zw3−∂2 Zw3+4νrh2w3 +(c0+cd−ca)−h∂Z(divxw)−1 ε2h∂Z(divXw)+h∇h·∂Zw+ ∇h·Z∂2 Zw−∂2 Zw3 =2νr ε2h2∂X1v2−∂X2v1+2νrh2∂x1v2−∂x2v1. 3.2. Asymptotic expansion Now we formally expand the unknowns: v(x,X,Z)=v0(x,X,Z)+εv1(x,X,Z)+ε2v2(x,X,Z)+ · · · ,(18) v3(x,X,Z)=v0 3(x,X,Z)+εv1 3(x,X,Z)+ε2v2 3(x,X,Z)+ · · · ,(19) p(x,X,Z)=1 ε2p0(x,X,Z)+1 εp1(x,X,Z)+p2(x,X,Z)+ · · · ,(20) w(x,X,Z)=w0(x,X,Z)+εw1(x,X,Z)+ε2w2(x,X,Z)+ · · · ,(21) w3(x,X,Z)=w0 3(x,X,Z)+εw1 3(x,X,Z)+ε2w2 3(x,X,Z)+ · · · .(22) From condition (13), we assume that the functions vi, vi 3,pi,wi, wi 3,i=1,2, . . . ,of(18)–(22), are T2-periodic functions in the variable X. I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 1921 The procedure is standard: we plug the above expansions into the rescaled equations (15)–(17) and collect the terms with equal powers of ε. For that purpose, we need to determine the asymptotic behavior of the terms involving function h. Taking into account (12)–(14), we deduce h2=ε2h2 1+2ε3h1h2+ε4h2 2∼O(ε2), h∇h=εh1∇Xh2+ε2h1∇xh1+ε2h2∇Xh2+ε3h2∇xh1∼O(ε), h1h=1 εh1∆Xh2+h2∆Xh2+ε2h1∆xh1+ε3h2∆xh1∼O1 ε, |∇h|2= |∇Xh2|2+2ε∇xh1∇Xh2+ε2|∇xh1|2∼O(1). 3.3. Main order term We start by substituting the expansions (18)–(22) into momentum and divergence equation (15)–(16). The leading order terms are given by 1 ε2: −(ν +νr)h2 1∆Xv0+h2 1∇Xp0=0, 1 ε2: −(ν +νr)h2 1∆Xv0 3=0, 1 ε:h1divXv0=0. (23) Note that (23)1,(23)3is, in fact, a Stokes system for (v0,p0)with respect to X. On the other hand, the third component v0 3 satisfies a simple Laplace equation (23)2, again with respect to X. Therefore, taking into account the boundary conditions with respect to X, we deduce ∇Xv0=0,∇Xp0=0,∇Xv0 3=0.(24) The main order term from the angular momentum equation (17) yields 1 ε2: −(ca+cd)h2 1∆Xw0−(c0+cd−ca)h2 1∇X(divXw0)=0, 1 ε2: −(ca+cd)h2 1∆Xw0 3=0. (25) Similarly as above, we conclude ∇Xw0=0,∇Xw0 3=0.(26) As we can see, main order terms led to a decoupled problem: (23) involves only the velocity and pressure, while (25) is satisfied only by the microrotation. Consequently, we established that the leading order terms v0, v0 3,w0, w0 3do not depend on the fast variable X. 3.4. Lower order terms We continue the computation and write the problems satisfied by the lower-order terms in the rescaled equations. In view of (24) and (26), from (15)–(16) we get 1 ε: −(ν +νr)h2 1∆Xv1+h2 1∇Xp1+(ν +νr)h1∆Xh2Z∂Zv0−h1∇Xh2Z∂Zp0=0, 1 ε: −(ν +νr)h2 1∆Xv1 3+(ν +νr)h1∆Xh2Z∂Zv0 3+h1∂Zp0=0, 1:h1divXv1− ∇Xh2·Z∂Zv0+∂Zv0 3=0. (27) Observe that there is no contribution of the terms involving microrotation field so we can proceed similarly as in [13]. We compute the mean value in Xof Eq. (27)3. Since h1,v0, v0 3do not depend on X, we find ∂v0 3=0. Using the boundary conditions for v0 3at the top and the bottom of the domain, we deduce that v0 3=0. Then the free-divergence condition written at order εgives h1divXv1= ∇Xh2·Z∂Zv0.(28) 1922 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 Taking the mean value in Xof Eq. (27)2, we conclude that ∂Zp0=0 and next ∇Xv1 3=0.(29) Taking the curlXoperator of the horizontal component in Eq. (27)1, we get h1curlX∆Xv1= ∇⊥ X∆Xh2·Z∂Zv0. Since h1and v0do not depend on X, h1curlXv1= ∇⊥ Xh2·Z∂Zv0.(30) We have ∆Xv1= ∇XdivXv1− ∇⊥ XcurlXv1and ∇X∇Xh2− ∇⊥ X∇⊥ Xh2=(∆Xh2)Id, thus using expressions (28) and (30), we obtain h2 1∆Xv1=h1∆Xh2Z∂Zv0.(31) Therefore, Eq. (27) can be written as ∇Xp1=0. The two next terms from the microrotation expansion are given by (see (17)): 1 ε:(ca+cd)∆Xw1+(c0+cd−ca)∇X(divXw1)=(c0+2cd)∆Xh2 h1 Z∂Zw0, 1 ε:(ca+cd)∆Xw1 3=(ca+cd)∆Xh2 h1 Z∂Zw0 3, (32) 1:(ca+cd)h2∆Xh2Z∂Zw0− |∇Xh2|2Z∂Zw0− |∇Xh2|2Z2∂2 Zw0−∂2 Zw0 +(c0+cd−ca)h2∆Xh2Z∂Zw0− |∇Xh2|2Z∂Zw0− |∇Xh2|2Z2∂2 Zw0 + ∇Xh2·∂Zw0 3+ ∇Xh2·Z∂2 Zw0 3+(ca+cd)−2h1h2∆Xw1+2h1∇Xh2Z∇X∂Zw1+h1∆Xh2Z∂Zw1 +(c0+cd−ca)−2h1h2∇X(divXw1)+h1∇Xh2Z∂Z(divXw1)+h1∆Xh2Z∂Zw1 +h1∇Xh2Z∇X∂Zw1−h1∇X∂Zw1 3−(ca+cd)h2 1∆Xw2−(c0+cd−ca)h2 1∇X(divXw2)=0, 1:(ca+cd)h2∆Xh2Z∂Zw0 3− |∇Xh2|2Z∂Zw0 3− |∇Xh2|2Z2∂2 Zw0−∂2 Zw0 3 +(c0+cd−ca)∇Xh2·Z∂2 Zw0−∂2 Zw0 3 +(ca+cd)−2h1h2∆Xw1 3+2h1∇Xh2·Z∇X∂Zw1 3+h1∆Xh2Z∂Zw1 3 −(c0+cd−ca)h1∂Z(divXw1)−(ca+cd)h2 1∆Xw2 3=0. (33) Let us prove that w0=0. For that purpose, we take the mean value with respect to Xin (33)1and carefully treat each term of this equation: (i) Terms involving w0, w0 3: since w0, w0 3do not depend on X, we have (ca+cd)T2h2∆Xh2Z∂Zw0− |∇Xh2|2Z∂Zw0− |∇Xh2|2Z2∂2 Zw0−∂2 Zw0dX +(c0+cd−ca)T2−h2∆Xh2Z∂Zw0− |∇Xh2|2Z∂Zw0− |∇Xh2|2Z2∂2 Zw0+ ∇Xh2·∂Zw0 3+ ∇Xh2·Z∂2 Zw0 3dX =(ca+cd)T2 h2∆Xh2dXZ∂Zw0−T2 |∇Xh2|2dXZ∂Zw0−T2 |∇Xh2|2dXZ2∂2 Zw0−∂2 Zw0 +(c0+cd−ca)−T2 h2∆Xh2dXZ∂Zw0−T2 |∇Xh2|2dXZ∂Zw0 −T2 |∇Xh2|2dXZ2∂2 Zw0+T2 ∇Xh2·∂Zw0 3dX+T2 ∇Xh2·Z∂2 Zw0 3dX =(ca+cd)T2−T2 |∇Xh2|2dXZ∂Zw0−T2 |∇Xh2|2dXZ∂Zw0−T2 |∇Xh2|2dXZ2∂2 Zw0−∂2 Zw0 +(c0+cd−ca)−T2 |∇Xh2|2dXZ∂Zw0−T2 |∇Xh2|2dXZ∂Zw0−T2 |∇Xh2|2dXZ2∂2 Zw0dX = −2(c0+2cd)MZ∂Zw0−(c0+2cd)MZ2∂2 Zw0−(ca+cd)∂2 Zw0. I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 1923 Here and in the sequel we introduce M=T2 |∇Xh2|2dX (34) as a coefficient depending on the considered rugosity profile. (ii) Using (32)1and the fact that w0does not depend on X, we obtain T2−2(ca+cd)h1h2∆Xw1−2(c0+cd−ca)h1h2∇X(divXw1)dX = −2T2 h1h2(ca+cd)∆Xw1+(c0+cd−ca)∇X(divXw1)dX = −2T2 (c0+2cd)h1h2 ∆Xh2 h1 Z∂Zw0dX = −2(c0+2cd)T2 h2∆Xh2dXZ∂Zw0 =2(c0+2cd)T2 |∇Xh2|2dXZ∂Zw0 =2(c0+2cd)MZ∂Zw0. (iii) The remaining terms involving w1, w1 3: employing again (32)1and the fact that h1and w0do not depend on X, we get (ca+cd)T22h1∇Xh2·Z∇X∂Zw1+h1∆Xh2Z∂Zw1dX +(c0+cd−ca)T2h1∇Xh2Z∂Z(divXw1)+h1∆Xh2Z∂Zw1+h1∇Xh2Z∇X∂Zw1−h1∇X∂Zw1 3dX =(ca+cd)T2−2h1h2·Z∂Z(∆Xw1)+h1h2Z∂Z(∆Xw1)dX +(c0+cd−ca)T2−h1h2Z∂Z∇X(divXw1)+h1h2Z∂Z(∆Xw1)−h1h2Z∂Z(∆Xw1)dX = − T2 h1h2Z∂Z(ca+cd)∆Xw1+(c0+cd−ca)∇X(divXw1)dX = − T2 h1h2Z∂Z(c0+2cd)∆Xh2 h1 Z∂Zw0dX = − T2 Z∂Z(c0+2cd)h2∆Xh2Z∂Zw0dX = −(c0+2cd)T2 h2∆Xh2dXZ∂Zw0−(c0+2cd)T2 h2∆Xh2dXZ2∂2 Zw0 =(c0+2cd)T2 |∇Xh2|2dXZ∂Zw0+(c0+2cd)T2 |∇Xh2|2dXZ2∂2 Zw0 =(c0+2cd)MZ∂Zw0+(c0+2cd)MZ2∂2 Zw0. (iv) Terms involving w2: integrating by parts, because the function h1depends only on xand the unknown w2is periodic in X, it follows −T2 (ca+cd)h2 1∆Xw2dX −T2 (c0+cd−ca)h2 1∇X(divXw2)dX =0. Adding all contributions (i)–(iv), we easily obtain −(ca+cd)∂2 Zw0+(c0+2cd)MZ∂Zw0=0. Combining the above equation with the corresponding boundary condition, namely w0|Z=0,1=0 finally gives w0=0. Proceeding analogously in (33)2, we derive the equation satisfied by w0 3: −(c0+2cd)∂2 Zw0 3+(ca+cd)MZ∂Zw0 3=0 1930 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 +Ω ε∆xh1+1 ε∆Xhε 2 h1+εhε 2 Z∂Zwε·φε−Ω ε2|∇xh1+1 ε∇Xhε 2|2 (h1+εhε 2)2Z∂Zwεφε −Ω ε2∇xh1+1 ε∇Xhε 2 h1+εhε 2 · ∇xwε∂Z(Zφε)+Ω ε2|∇xh1+1 ε∇Xhε 2|2 (h1+εhε 2)2∂wε·∂Z(Z2φε) −Ω ∇xh1+1 ε∇Xhε 2 (h1+εhε 2)2·ε∂Zw3,εφε+Ω 1 (h1+εhε 2)2∇xw3,ε∂Zφε−Ω ∇xh1+1 ε∇Xhε 2 (h1+εhε 2)2ε∂Zw3,ε∂Z(Zφε) =2νrΩ 1 h1+εhε 2 ε−∂Zv2,ε i+∂Zv1,ε jφε+2νrΩ ε2rotxv3,εφε. Passing to the limit, we get (ca+cd)ΩT2 ∇Xw1· ∇Xφ+(ca+cd)ΩT2 ∆Xh2 h1 Z∂Zw0φ ×(c0+cd−ca)ΩT2 divXw1divXφ+(c0+cd−ca)ΩT2 ∆Xh2 h1 Z∂Zw0φ=0, for any φ∈D(Ω;C1(T2)). This is equivalent to (61). We proceed analogously with Eq. (53) in order to deduce the relation (62). Finally, from Lemma 4 we can easily get the additional result for the pressure: Lemma 6. The pressure is such that limε→0ΩεpεdivX(φ) =0for any φ∈L2(Ω;H1(T2)). 4.3.2. Passing to the limit Now we are in position to pass to the limit in the rescaled equations (49)–(53). Following same arguments as in [13, Sec. 5.2 and 5.3], from divergence and momentum equation it is straightforward to obtain the weak formulations corresponding to (57)–(59). It remains to verify Eq. (60). Let us start with the equation for w3,ε to confirm that w1 3=0. We employ εφ(x,Z),φ∈D(Ω)as a test function in (53). Using the identity 1 h1h−1 h2|∇h|2=div 1 h∇h, we have (ca+cd)Ω ε∇xw3,ε · ∇xφε−(ca+cd)Ω ε h1+εhε 2∇xh1+1 ε∇Xh2· ∇xw3,ε∂Z(Zφ) +(ca+cd)Ω ε|∇xh1+1 ε∇Xhε 2|2 (h1+εhε 2)2∂Zw3,ε ·∂Z(Z2φε)+(ca+cd)Ω 1 (h1+εhε 2)2 1 ε∂Zw3,ε ·∂Zφε ×4νrΩ εw3,εφε+(c0+cd−ca)Ω 1 h1+εhε 2 divxwε∂Zφε−(c0+cd−ca)Ω ε∇xh1+1 ε∇Xhε 2 h1+εhε 2 ·wε∂Z(φε) −(c0+cd−ca)Ω ε∇xh1+1 ε∇Xhε 2 h1+εhε 2 ·1 ε∂Zwε∂Z(Zφε)+(c0+cd−ca)Ω 1 (h1+εhε 2)2 1 ε∂Zw3,ε∂Zφε =2νrΩ ε∂x1v2,ε −∂x2v1,εφε. Passing to the limit and taking into account that w1 3=w1 3(x,Z), we obtain −(ca+cd)ΩT2 1 h1 ∇Xh2∇Xw2 3∂Z(Zφ) +(ca+cd)ΩT2 |∇Xh2|2 h2 1 ∂Zw1 3∂Z(Z2φ) +(ca+cd)ΩT2 1 h2 1 ∂Zw1 3∂Zφ+(c0+cd−ca)ΩT2 1 h2 1 ∂Zw1 3∂Zφ=0. Using relation (62) from Lemma 5, it can be easily verified that the above relation is, in fact, the energy formulation corresponding to −(ca+cd)T2 h2∆Xh2dX1 h2 1 Z∂Z(Z∂Zw1 3)−(ca+cd)T2 |∇Xh2|2dX1 h2 1 Z2∂Z(∂Zw1 3) −(c0+2cd)T2 1 h2 1 ∂2 Zw1=0. I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 1931 Integrating by parts gives −(c0+2cd)T2 1 h2 1 ∂2 Zw1+(ca+cd)T2 |∇Xh2|2dX1 h2 1 Z2∂Z(∂Zw1 3)=0. Observe that this corresponds to Eq. (47) obtained in a formal way. In view of the zero boundary condition for microrotation, we conclude w1 3=0. Let us proceed with equation for wε. Analogously, we multiply Eq. (52) by εφ(x,Z),φ∈D(Ω)to obtain (ca+cd)Ω ε∇xwε· ∇xφε−(ca+cd)Ω ε h1+εhε 2∇xh1+1 ε∇Xh2· ∇xwε∂Z(Zφ) +(ca+cd)Ω ε|∇xh1+1 ε∇Xhε 2|2 (h1+εhε 2)2∂Zwε·∂Z(Z2φε)+(ca+cd)Ω 1 (h1+εhε 2)2 1 ε∂Zwε·∂Zφε ×4νrΩ εwεφε+(c0+cd−ca)Ω divxwε(εdivxφε)−(c0+cd−ca)Ω ε∇xh1+1 ε∇Xhε 2 h1+εhε 2 ·divxwε∂Z(Zφε) +(c0+cd−ca)Ω ε|∇x+1 ε∇Xhε 2|2 (h1+εhε 2)2∂Zwε·∂Z(Z2φε)+(c0+cd−ca)Ω ε∇x+1 ε∇Xh2 (h1+εhε 2)2 1 ε∂Zw3,εφ +(c0+cd−ca)Ω 1 (h1+εhε 2)2 1 ε∂Zw3,ε∇xφ−(c0+cd−ca)Ω ∇x+1 ε∇Xh2 (h1+εhε 2)2 1 ε∂Zw3,ε∂Z(Zφ) =2νrΩ 1 h1+εhε 2−∂Zv2,ε i+∂Zv1,ε jφε+2νrΩ εrotxv3,εφε. Passing to the limit and taking into account that w1=w1(x,Z)and w1 3=0 give −(ca+cd)ΩT2 1 h1 ∇Xh2· ∇Xw2∂Z(Zφ) +(ca+cd)ΩT2 |∇Xh2|2 (h1)2∂Zw1∂Z(Z2φ) +(ca+cd)ΩT2 1 (h1)2∂Zw1·∂Zφ−(c0+cd−da)ΩT2 1 h1 ∇Xh2·divXw2∂Z(Zφ) +(c0+cd−ca)ΩT2 |∇Xh2|2 (h1)2∂Zw1∂Z(Z2φ) =2νrΩT2 1 h1 (−∂Zv0 2i+∂Zv0 1j)φ. Using (61) from Lemma 5, it can be easily verified that the latter is the energy formulation corresponding to −(c0+2cd)T2 h2∆Xh21 h2 1 Z∂Z(Z∂Zw1)−(ca+cd)T2 |∇Xh2|2dX1 h2 1 Z2∂Z(∂Zw1) −(c0+cd−ca)T2 |∇Xh2|21 h2 1 Z2∂Z(∂Zw1)−(ca+cd) h2 1 ∂2 Zw1=2νr h1 (−∂Zv0 2i+∂Zv0 1j). Finally, after integrating by parts, we get −(ca+cd)∂2 Zw1+(c0+2cd)T2 |∇Xh2|2dXZ∂Zw1=2νrh1(−∂Zv0 2i+∂Zv0 1j) which corresponds to (60). Acknowledgments The first author of this work has been supported by the Croatian Science Foundation (scientific project 3955: Mathematical modeling and numerical simulations of processes in thin or porous domains). The second author has been supported by the projects MTM2011-24457 of the Ministerio de Economía y Competitividad and FQM309 of the Junta de Andalucía. References [1] O. Reynolds, On the theory of lubrication and its applications to Mr. Beauchamp Tower’s experiments, including an experimental determination of the viscosity of olive oil, Philos. Trans. R. Soc. Lond. 177 (1886) 157–234. [2] G.H. Wannier, A contribution to the hydrodynamics of lubrication, Quart. Appl. Math. 8 (1950) 1–32. [3] H.G. Elrod, A derivation of the basic equations for hydrodynamics lubrication with a fluid having constant properties, Quart. Appl. Math. 17 (1960) 349–359. [4] G. Bayada, M. Chambat, The transition between the Stokes equations and the Reynolds equation: a mathematical proof, Appl. Math. Optim. 14 (1986) 73–93. [5] G.J. Jonhnston, R. Wayte, H.A. Spikes, The measurement and study of very thin lubricant films in concentrated contacts, Tribol. Trans. 34 (1991) 187–194. [6] J.B. Luo, P. Huang, S.Z. Wen, Thin film lubrication part I: study on the transition between EHL and thin film lubrication using relative optical interference intensity technique, Wear 194 (1996) 107–115. 1932 I. Pažanin, F.J. Suárez-Grau / Computers and Mathematics with Applications 68 (2014) 1915–1932 [7] J.B. Luo, P. Huang, S.Z. Wen, L. Lawrence, Characteristics of fluid lubricant films at nano-scale, J. Tribol. 121 (1999) 872–878. [8] A.C. Eringen, Theory of micropolar fluids, J. Math. Mech. 16 (1966) 1–16. [9] G. Lukaszewicz, Micropolar Fluids: Theory and Applications, Birkhäuser, Boston, 1999. [10] A. Dyson, Hydrodynamic lubrication of rough surface—a review work, in: Proceedings of the 4th Leeds–Lyon Symposium on Surfaces Roughness on Lubrication, 1977, pp. 61–69. [11] G. Bayada, M. Chambat, New models in the theory of the hydrodynamic lubrication of rough surfaces, J. Tribol. 110 (1988) 402–407. [12] N. Benhaboucha, M. Chambat, I. Ciuperca, Asymptotic behaviour of pressure and stresses in a thin film flow with a rough boundary, Quart. Appl. Math. 63 (2005) 369–400. [13] D. Bresch, C. Choquet, L. Chupin, T. Colin, M. Gisclon, Roughness-induced effect at main order on the Reynolds approximation, SIAM Multiscale Model. Simul. 8 (2010) 997–1017. [14] L. Chupin, S. Martin, Rigorous derivation of the thin film approximation with roughness-induced correctors, SIAM J. Math. Anal. 44 (2012) 3041–3070. [15] J.-L. Ligier, Lubrification des Paliers Moteurs, Technip, 1997. [16] J. Casado-Diáz, M. Luna-Laynez, F.J. Suárez-Grau, Asymptotic behavior of the Navier–Stokes system in a thin domain with Navier condition on a slightly rough boundary, SIAM J. Math. Anal. 45 (2013) 1641–1674. [17] I. Pažanin, F.J. Suárez-Grau, Effects of rough boundary on the heat transfer in a thin-film flow, C. R. Mec. 341 (2013) 646–652. [18] M.C. Pereira, Parabolic problems in highly oscillating thin domains, Ann. Mat. Pura Appl. (2014) http://dx.doi.org/10.1007/s10231-014-0421-7. [19] C. Singh, P. Sinha, The three-dimensional Reynolds’ equation for micropolar fluid lubricated bearings, Wear 76 (1982) 199–209. [20] G. Bayada, G. Lukaszewicz, On micropolar fluids in the theory of lubrication. Rigorous derivation of an analogue of the Reynolds equation, Internat. J. Engrg. Sci. 34 (1996) 1477–1490. [21] D. Dupuy, G. Panasenko, R. Stavre, Asymptotic solution for a micropolar flow in a curvilinear channel, ZAMM Z. Angew. Math. Mech. 88 (2008) 793–807. [22] E. Marušić-Paloka, I. Pažanin, S. Marušić, An effective model for the lubrication with micropolar fluid, Mech. Res. Comm. 52 (2013) 69–73. [23] M. Boukrouche, L. Paoli, Asymptotic analysis of a micropolar fluid flow in a thin domain with a free and rough boundary, SIAM J. Math. Anal. 44 (2012) 1211–1256. [24] G. Bayada, N. Benhaboucha, M. Chambat, New models in micropolar fluid and their applications to lubrication, Math. Models Methods Appl. Sci. 15 (2005) 343–374. [25] C. Conca, F. Murat, O. Pironneau, The Stokes and Navier–Stokes equations with boundary conditions involving the pressure, Jpn. J. Math. 20 (1994) 263–318. [26] G. Nguetseng, A general convergence result for a functional related to the theory of homogenization, SIAM J. Math. Anal. 20 (1989) 608–623. [27] G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23 (1992) 1482–1518. [28] S. Marušić, E. Marušić-Paloka, Two-scale convergence for thin domains and its applications to some lower-dimensional models in fluid mechanics, Asymptot. Anal. 23 (2000) 23–58.