# Analysis of the Reynolds Equation for Lubrication in Case of Pressure-Dependent Viscosity

Authors: Eduard Marušić-Paloka; Sanja Marušić
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

We study the Reynolds equation, describing the  ow of a lubricant, in case of pressure-dependent viscosity. First we prove the existence and uniqueness of the solution. Then, we study the asymptotic behavior of the solution in case of periodic roughness via homogenization method. Some interesting nonlocal effects appear due to the nonlinearity.


---

## Body

## 1. Introduction

The Reynolds equations [1] describe the flow of a thin film of lubricant separating two rigid surfaces in relative motion. Controlling the flow of lubricant is an important engineering issue since inappropriate lubrication would increase the friction and wear, finally resulting in degrading the performance of the device. In his classical paper from 1848 [2] Stokes predicted that the viscosity of the fluid can depend on the pressure. Those effects for various liquids have been measured in many engineering papers, starting from the beginning of the 20th century (see, e.g., [3]). That effect is usually neglected as it becomes important only in case of high pressure. Most fluid-lubricated bearings operate with high pressure and in such a flow regime the dependence of the viscosity on the pressure becomes important. According to Szeri [4] the idea of pressure-dependent viscosity was introduced in lubrication theory by Gatcombe in 1945 [5]. Several models have been used to describe that relation since. The most popular is probably the exponential law
(1)μ=μ0exp(αp)
usually called the Barus formula [6]. Here μ0 and α are the constants depending on the lubricant. The formula seems to be reasonable for mineral oil, unless the pressure is very high (larger then 0.5 MPa). The coefficient α typically ranges between 1 and 10−8. The lower end of the range corresponds to paraffinic and the upper end corresponds to the naphthenic oils (see Jones [7]). That formula is still frequently used by engineers. The simplest viscosity-pressure relation is given by the the power law
(2)μ(p)=α|p|γ.
In case of the two above-mentioned laws explicit solutions of the equations of motion, for some particular situations like unidirectional and plane-parallel flows, were found in [8]. Discussion on other possibilities for the viscosity-pressure formula and some historical remarks on the subject can be found in the same paper. Several engineering papers can be found discussing other possible laws and their consistency. We mention, for instance, [9, 10].We do not make any assumption on the particular form of the functionp↦μ(p). Some technical assumptions, like smoothness, will be needed for the proofs.We study the stationary version of the Reynolds equation. Unless the velocity of relative motion is time dependent, steady approximation is reasonable in most cases (see, e.g., [4, Chapter 2.2]).Our first goal is to prove that the problem is well posed. Secondly, we investigate the asymptotic behavior of the solution in case of periodically distributed asperities. Using the homogenization approach, we find the macroscopic Reynolds pressure. Interesting nonlocal effects appear due to the nonlinearity caused by the pressure-dependent viscosity.


## 2. Position of the Problem

The fluid domain is bounded by two rigid surfaces. The simplified mathematical model can be written in the following form. Let𝒪⊂Rn be a bounded domain and let h:𝒪¯→R be a bounded, strictly positive smooth function such that
(3)0<h0≤h(x)≤h1.
Function h describes the shape of the slide. By ε≪1 we denote a very small parameter representing the domain thickness. Using the shape function h we define the fluid domain (Figure 1) by
(4)Ωε={X=(x,xn+1)∈Rn+1;mmXx=(x1,…,xn)∈𝒪,0<xn+1<εh(x)Rn+1vdxn+1}.
We then consider the stationary flow through a domain Ωε. We want to describe the situation with a lower-dimensional model. The velocity of the relative motion of two surfaces is the constant vector denoted by V=v1e1+⋯+vnen. The unknowns in the model are u (the velocity) and p (the pressure). We recall that the stationary motion of the incompressible viscous laminar flow is governed by the stationary Navier-Stokes equations. Thus we write the following system:
(5)-div[μ(p)Du]+(u∇)u+∇p=0,divu=0inΩε,u=Vforxn+1=0,u=0forxn+1=εh(x),
where Du=(1/2)(∇u+∇ut) is the symmetric part of the velocity gradient. It is important to notice that in such system the pressure is not defined only up to a constant, as in the classical Navier-Stokes system with constant viscosity. Under certain technical assumptions, if the given data are not too large, the existence of the solution for such system was discussed in [11, 12]. Neglecting the effects of inertia, we get the Stokes system with pressure-dependent viscosity
(6)div[μ(p)Du]+∇p=0,divu=0inΩε,
studied in [13].Figure 1If the thickness of the domain is small, the solution can be fairly approximated by the solution of the Reynolds equations [4, 14]
(7)v(x,xn+1)m=112μ(p)xn+1(εh-xn+1)∇p(x)+(1-xn+1εh)V,(8)div(∫0εh‍vdxn+1)=0.
Indeed, if we derive a formal asymptotic expansion of the solution to the system (5) in powers of ε, then the solution of the Reynolds equation (7) makes the first term of the expansion (see, e.g., [4]). Here and in the sequel the differential operatorsdivand ∇ are taken only with respect to x variable; that is,
(9)divb=∂b1∂x1+⋯+∂bn∂xn,∇ϕ=∂ϕ∂x1e1+⋯+∂ϕ∂xnen.
It leads to an elliptic equation of the form
(10)div(h3μ(p)∇p)=6V·∇hin𝒪(11)p=qon∂𝒪.
The goal of this paper is to study that equation.We assume that the functionp↦μ(p) is of class C1(R) and μ>0 for any value of p. In real life the viscosity increases with pressure, but such an assumption is not necessary for our study.


## 3. Existence of the Solution

### 3.1. Transformed Equation

Equation (10) is a quasilinear elliptic PDE, but it can be linearized by simple trick. To do so we rewrite the equation using the function
(12)Mσ(p)=∫σp‍dtμ(t).
We choose σ≤q. Function Mσ is strictly increasing, since Mσ′(p)=(1/μ(p))>0 and thus it is bijective. Furthermore Mσ(p) has the same sign as p-σ; that is, for p>σ we have Mσ(p)>0, for p<σ obviously Mσ(p)<0, and finally Mσ(p)=0, if and only if p=σ.We introduce the new unknown function(13)w(x)=Mσ(p(x))=∫σp(x)‍dtμ(t).
At this point we assume that the integral ∫0-∞‍(dt/μ(t)) is divergent, that is,
(14)∫-∞0‍dtμ(t)=+∞.
As a consequence
(15)limσ→-∞Mσ(p)=+∞,∀p∈R,
as well as
(16)limp→-∞Mσ(p)=-∞,∀σ∈R.
Deriving (13) we obtain
(17)1μ(p)∇p=∇w
and the problem can be written as
(18)div(h3∇w)=6V·∇hin𝒪,(19)w=Mσ(q)on∂𝒪.
That is a linear elliptic equation for w and it has a unique solution. To get the existence and uniqueness of the solution we quote Theorem 8.34 from classical book of Gilbarg and Trudinger [15]. For simplicity, here and in the sequel, we assume that q and, consequently, Mσ(q) are defined on whole 𝒪¯. We combine that with the maximum principle from the appendix, and it gives the following.Theorem 1.
Under the assumption that the boundary ∂𝒪 is of class C1,α and that h∈Cα(𝒪¯), q∈C1,α(𝒪¯) the problem (18), (19) has a unique solution
(20)w∈C1,α(𝒪¯).
Furthermore
(21)w(x)≤M(q¯)+𝒵,
where
(22)q¯=|q|L∞(∂ω)
and 𝒵=0 if V·∇h<0. Otherwise
(23)𝒵={6|V|[∫01dth(t)2-h03h15]+(h0h1)3∫q0q1dsμ(s)ifn=13(8/5)(32)(28/5)(2π)(1/4)h03×d(6|V|h1|𝒪|1/5+h13|∇Mσ(q)|L5(𝒪))|𝒪|(7/5)ifn=2,
withd=diam𝒪.Proof.
The existence follows directly from Theorem 8.34. from Gilbarg and Trudinger [15]. If V·∇h<0, then (21) follows directly from the weak maximum principle (see, e.g., [15]). In case n=1 the problem can be solved by quadratures and the solution given by (25) can be easily estimated to get (23). In the remaining case n=2 (21) follows from the special variant of the maximum principle proved in the appendix.Remark 2.
In casen=1 (18) is an ODE (we take 𝒪=]0,1[  and q0>q1, without losing generality)
(24)(h3w′)′=6Vh′in]0,1[,w(i)=Mσ(qi)fori=0,1,
and it can be solved by quadratures
(25)w(x)=Mσ(q0)m+∫0x‍1h(t)3[6V(h(t)-∫01‍(dr/h(r)2)∫01‍(ds/h(s)3))mmmmmmmm+Mσ(q1)-Mσ(q0)∫01‍(ds/h(s)3)]dt.


### 3.2. Back to the Original Equation

Now, our goal is not to find the auxiliary functionw but to find the pressure p. Since we have introduced w as
(26)w=Mσ(p),
we should have p=Mσ-1(w). In order to do so we have to make sure that w(x)∈ImMσ for any x∈𝒪¯. Since Mσ is strictly increasing and we have assumed that (14) holds, if we define
(27)Mσ+=lims→+∞Mσ(s)=∫σ+∞‍dsμ(s),
due to (16) we obviously have for any p∈R(28)-∞=lims→-∞Mσ(s)≤Mσ(p)≤Mσ+.
Thus
(29)ImMσ=]-∞,Mσ+[.
So, to fulfill the condition w(x)∈ImMσ we need to have
(30)w(x)≤Mσ+,∀x∈𝒪¯.That condition is not necessarily fulfilled.In view of (21) that condition reduces to
(31)𝒵≤∫q¯+∞‍dtμ(t),
where 𝒵=𝒵(V,h,𝒪,μ,q) is defined by (23).We have proved the following theorem.Theorem 3.
Suppose that the conditions of Theorem1 hold, and that in addition (31) is fulfilled. Then p=Mσ-1(w) is the unique solution of (10) and (11).Remark 4.
It is important to notice that even thoughw  does depend on σ, the effective pressure p  does not. For the purpose of this remark, we denote Hσ(w)=Mσ-1(w)   to stress the dependence on the parameter σ  which is of interest here. We start by
(32)p=Hσ(w).
It is obvious from the definition of Mσ  that
(33)∂Mσ∂p(p)=1μ(p),∂Mσ∂σ(p)=-1μ(σ).
As Mσ(Hσ(w))=w, deriving with respect to σ we arrive at
(34)∂Hσ∂σ(w)=-(∂Mσ/∂σ)(p)(∂Mσ/∂p)(p)=μ(p)μ(σ)=μ(Hσ(w))μ(σ).
Deriving (13) we get (∂w/∂σ)=-1/μ(σ). Using the rule for deriving the inverse function, we have
(35)∂Hσ∂w(w)=μ(p)=μ(Hσ(w)).
Thus
(36)∂p∂σ=∂Hσ∂w(w)∂w∂σ(σ)+∂Hσ∂σ(w)=-μ(Hσ(w))μ(σ)+μ(Hσ(w))μ(σ)=0.


## 3.1. Transformed Equation

Equation (10) is a quasilinear elliptic PDE, but it can be linearized by simple trick. To do so we rewrite the equation using the function
(12)Mσ(p)=∫σp‍dtμ(t).
We choose σ≤q. Function Mσ is strictly increasing, since Mσ′(p)=(1/μ(p))>0 and thus it is bijective. Furthermore Mσ(p) has the same sign as p-σ; that is, for p>σ we have Mσ(p)>0, for p<σ obviously Mσ(p)<0, and finally Mσ(p)=0, if and only if p=σ.We introduce the new unknown function(13)w(x)=Mσ(p(x))=∫σp(x)‍dtμ(t).
At this point we assume that the integral ∫0-∞‍(dt/μ(t)) is divergent, that is,
(14)∫-∞0‍dtμ(t)=+∞.
As a consequence
(15)limσ→-∞Mσ(p)=+∞,∀p∈R,
as well as
(16)limp→-∞Mσ(p)=-∞,∀σ∈R.
Deriving (13) we obtain
(17)1μ(p)∇p=∇w
and the problem can be written as
(18)div(h3∇w)=6V·∇hin𝒪,(19)w=Mσ(q)on∂𝒪.
That is a linear elliptic equation for w and it has a unique solution. To get the existence and uniqueness of the solution we quote Theorem 8.34 from classical book of Gilbarg and Trudinger [15]. For simplicity, here and in the sequel, we assume that q and, consequently, Mσ(q) are defined on whole 𝒪¯. We combine that with the maximum principle from the appendix, and it gives the following.Theorem 1.
Under the assumption that the boundary ∂𝒪 is of class C1,α and that h∈Cα(𝒪¯), q∈C1,α(𝒪¯) the problem (18), (19) has a unique solution
(20)w∈C1,α(𝒪¯).
Furthermore
(21)w(x)≤M(q¯)+𝒵,
where
(22)q¯=|q|L∞(∂ω)
and 𝒵=0 if V·∇h<0. Otherwise
(23)𝒵={6|V|[∫01dth(t)2-h03h15]+(h0h1)3∫q0q1dsμ(s)ifn=13(8/5)(32)(28/5)(2π)(1/4)h03×d(6|V|h1|𝒪|1/5+h13|∇Mσ(q)|L5(𝒪))|𝒪|(7/5)ifn=2,
withd=diam𝒪.Proof.
The existence follows directly from Theorem 8.34. from Gilbarg and Trudinger [15]. If V·∇h<0, then (21) follows directly from the weak maximum principle (see, e.g., [15]). In case n=1 the problem can be solved by quadratures and the solution given by (25) can be easily estimated to get (23). In the remaining case n=2 (21) follows from the special variant of the maximum principle proved in the appendix.Remark 2.
In casen=1 (18) is an ODE (we take 𝒪=]0,1[  and q0>q1, without losing generality)
(24)(h3w′)′=6Vh′in]0,1[,w(i)=Mσ(qi)fori=0,1,
and it can be solved by quadratures
(25)w(x)=Mσ(q0)m+∫0x‍1h(t)3[6V(h(t)-∫01‍(dr/h(r)2)∫01‍(ds/h(s)3))mmmmmmmm+Mσ(q1)-Mσ(q0)∫01‍(ds/h(s)3)]dt.


## 3.2. Back to the Original Equation

Now, our goal is not to find the auxiliary functionw but to find the pressure p. Since we have introduced w as
(26)w=Mσ(p),
we should have p=Mσ-1(w). In order to do so we have to make sure that w(x)∈ImMσ for any x∈𝒪¯. Since Mσ is strictly increasing and we have assumed that (14) holds, if we define
(27)Mσ+=lims→+∞Mσ(s)=∫σ+∞‍dsμ(s),
due to (16) we obviously have for any p∈R(28)-∞=lims→-∞Mσ(s)≤Mσ(p)≤Mσ+.
Thus
(29)ImMσ=]-∞,Mσ+[.
So, to fulfill the condition w(x)∈ImMσ we need to have
(30)w(x)≤Mσ+,∀x∈𝒪¯.That condition is not necessarily fulfilled.In view of (21) that condition reduces to
(31)𝒵≤∫q¯+∞‍dtμ(t),
where 𝒵=𝒵(V,h,𝒪,μ,q) is defined by (23).We have proved the following theorem.Theorem 3.
Suppose that the conditions of Theorem1 hold, and that in addition (31) is fulfilled. Then p=Mσ-1(w) is the unique solution of (10) and (11).Remark 4.
It is important to notice that even thoughw  does depend on σ, the effective pressure p  does not. For the purpose of this remark, we denote Hσ(w)=Mσ-1(w)   to stress the dependence on the parameter σ  which is of interest here. We start by
(32)p=Hσ(w).
It is obvious from the definition of Mσ  that
(33)∂Mσ∂p(p)=1μ(p),∂Mσ∂σ(p)=-1μ(σ).
As Mσ(Hσ(w))=w, deriving with respect to σ we arrive at
(34)∂Hσ∂σ(w)=-(∂Mσ/∂σ)(p)(∂Mσ/∂p)(p)=μ(p)μ(σ)=μ(Hσ(w))μ(σ).
Deriving (13) we get (∂w/∂σ)=-1/μ(σ). Using the rule for deriving the inverse function, we have
(35)∂Hσ∂w(w)=μ(p)=μ(Hσ(w)).
Thus
(36)∂p∂σ=∂Hσ∂w(w)∂w∂σ(σ)+∂Hσ∂σ(w)=-μ(Hσ(w))μ(σ)+μ(Hσ(w))μ(σ)=0.


## 4. Homogenization

In this section we want to study the effects of rugosities of surfaces on lubrication process. The idea of finding the macroscopic effects of roughness on lubrication process, via homogenization, is quite old and well studied. Case of constant viscosity for incompressible and compressible flows as well as non-Newtonian, deformation dependent, viscosities were investigated. The subject was treated by several authors and we here mention [16–18]. The case of pressure-dependent viscosity brings some new interesting nonlocal effects.We assume that the functionh, describing the form of the fluid domain, is periodic with small period 1/m, with m∈N. To stress that dependence we denote it by hm. More precisely, we denote by Y=]0,1[n,n=1,2, the period. We further assume that h:Rn→[d0,+∞[, d0>0 is periodic with period Y and smooth. Then we take hm of the form
(37)hm(x)=h(mx).
Thus, the function h describes the form of periodically distributed rugosities.To emphasize that the relative velocity of bearing surfacesV is large, we assume that it also depends on m, the same parameter that is taken for description of rugosities. In that case our equation reads
(38)div(hm3μ(pm)∇pm)=6Vm·∇hmin𝒪.


### 4.1. One-Dimensional Case

Ifn=1, the above problem is posed on an interval 𝒪=]0,1[. With an appropriate boundary condition
(39)pm(0)=q0,pm(1)=q1.
It forms a boundary value problem for nonlinear ODE:
(40)(hm3μ(pm)pm′)′=6Vmhm′in𝒪.
To study the asymptotic behavior of the solution with respect to m we linearize the problem using the transformation wm=Mσ(pm). To simplify, in this section we choose σ=q0 and, dropping the index σ in Mσ and Mσ+, we denote
(41)M(p)=∫q0p‍dtμ(t),M+=∫q0+∞‍dtμ(t).Theorem 5.
Let
(42)
χ
0
(
y
)
=
6
(
∫
0
1
‍
d
s
h
(
s
)
3
)
-
1
m
×
[
〈
1
h
3
〉

∫
0
y
‍
d
v
h
(
v
)
2
-

〈
1
h
2
〉

∫
0
y
‍
d
v
h
(
v
)
3
]
,
(43)
〈
·
〉
=
∫
0
1
‍
·
d
y
,
and let
H
=
M
-
1. Suppose that there exists a limit(44)V=limm→∞Vmmand that, form  large enough andM+  defined in (27), the following condition holds:(45)b3m(x)b3m(1)M(q1)+6Vmb3m(1)b2m(x)-b3m(x)b2m(1)b3m(1)≤M+,with(46)bαm(x)=∫0x‍dthm(t)α,α=2,3.Then(47)pm⇀p0=∫01‍H(xM(q1)+Vχ0(y))dymmmmmmmweak*inL∞(0,1).Proof.
Equation (40), with boundary conditions (39), can be solved by reduction to quadratures, after the substitution
(48)wm=M(pm),
with Mσ strictly increasing function defined by (12). The problem for wm now reads
(49)(hm3wm′)′=6Vmhm′in]0,1[,wm(0)=M(q0)=0,wm(1)=M(q1).
It is easy to see that (49) has a unique solution given by (25). Since M(q0)=0, (25) now reduces to
(50)wm=b3m(x)b3m(1)Mσ(q1)+6Vmb3m(1)b2m(x)-b3m(x)b2m(1)b3m(1),
where, for α=2,3,
(51)bαm(x)=∫0x‍dthm(t)α=1m∫0mx‍dsh(s)α⟶x∫01‍dsh(s)α,mmmmmmmmmasm⟶+∞.
Now pm, the solution to the problem (40), exists if (45) is fulfilled. The second term in (50) thus obviously tends to xM(q1), as m→+∞. The last term is more interesting. The denominator tends to
(52)〈1h3〉=∫01‍dsh(s)3.
As for its numerator, we have
(53)Vm[b3m(1)b2m(x)-b3m(x)b2m(1)]mm=Vm[∫01‍∫0x‍(1hm(t)3hm(s)2-1hm(t)2h(s)m3)dsdt]mm=Vmm(〈1h3〉∫0mx‍dvh(v)2-〈1h2〉∫0mx‍dvh(v)3).
Suppose that(54)V=limm→+∞Vmm
and denote
(55)G(y)=〈1h3〉∫0y‍dvh(v)2-〈1h2〉∫0y‍dvh(v)3.
Obviously the function G is periodic with period 1 so that, due to the standard periodicity lemma (see, e.g., [19]), as m→+∞,
(56)G(mx)⇀〈G〉=∫01‍G(y)dy,weak*inL∞(0,1).
By direct computation(57)〈G〉=〈yh3〉〈1h2〉-〈1h3〉〈yh2〉=∫∫01s-th(s)3h(t)2dsdt.
Thus
(58)Vm[b3m(1)b2m(x)-b3m(x)b2m(1)]⇀V〈G〉.
Now, denoting(59)χ0(y)=G(y)b3(1)
we have
(60)wm≈xM(q1)+Vχ0(mx),
and thus
(61)wm⇀w0=xM(q1)+V〈χ0〉weak*inL∞(0,1).
However, we are not interested in convergence of the auxiliary functionwm but in the convergence of the pressure pm. Since (45) is assumed to be true, we can define pm=H(wm), where H=M-1, and we have
(62)pm⇀p0=∫01‍H(xM(q1)+Vχ0(y))dymmmmmweak*inL∞(0,1).Deriving the expression on the right-hand side, we obtain the effective pressure drop in the form(63)p0′(x)=M(q1)∫01‍μ(xM(q1)+Vχ0(y))dy(64)=∫q0q1‍dτμ(τ)∫01‍μ(xM(q1)+Vχ0(y))dy=∫q0q1‍dτμ(τ)∫01‍μ(x∫q0q1‍dsμ(s)+Vχ0(y))dy.
As we can see, the pressure drop is not constant, as for the Newtonian flow. The interesting effect appears if V≠0 because, in that case, the expressions for the pressure and for the pressure drop are nonlocal due to the integral with respect to y. That phenomenon is entirely due to the fact that the viscosity is depending on the pressure.


### 4.2. Two-Dimensional Case

We suppose here that the functionhm is constructed from positive, smooth, Y-periodic function h:R2→[d0,+∞[, d0>0, Y=]0,1[2 in the same way as before, that is, by taking
(65)hm(x)=h(mx).
We have seen in the previous section that interesting effects happen only if we assume that
(66)Vm=mV.
In that case our equation reads
(67)div(hm3μ(pm)∇pm)=6Vm·∇hmin𝒪.
As we did in the existence analysis and in the previous section, we linearize the equation by substitution
(68)wm=Mσ(pm),
where the function Mσ is defined by (12). Now wm satisfies
(69)div(h3M∇wm)=6Vm·∇hmin𝒪.We postulate the asymptotic expansion in the form(70)wm(x)≈w0(x,y)+1mw1(x,y)m+1m2w2(x,y)+⋯,y=mx.
All functions are assumed to be Y-periodic in y variable.Plugging that in (69) and collecting, formally, terms with equal powers of m, we getm2: divy(h3∇yw0)=6V·∇yhm: divy(h3∇yw1)+divy(h3∇xw0)+h3divx∇yw0=01:divy(h3∇yw2)+divy(h3∇xw1)+h3divx∇yw1+h3Δxw0=0.Denoting(71)V=|V|,V0=VV
we have
(72)w0(x,y)=v0(x)+Vχ0(y),divy(h3∇yχ0)=6V0·∇yhinY,w1(x,y)=∑k=12‍χk(y)∂v0∂xk(x)+v1(x),divy(h3∇y(χk+yk))=0,∑k,ℓ=12‍akℓ∂2v0∂xk∂xℓ=0in𝒪,akℓ=∫Y‍h3∂∂yℓ(χk+yk)dy.Remark 6.
The same computation can be done in one-dimensional case and it gives(73)(h3χ0′)′=6h′⇒χ0=6∫0y‍dth(t)2+C0∫0y‍dth(t)3+C1,(74)v0=xM(q1).
Constants C0,C1 are chosen in a way that boundary conditions χ0(0)=χ0(1)=0  are met, and it follows that
(75)C1=0,C0=-6(∫01‍dth(t)3)-1∫01‍dth(t)2.
Then
(76)wm≈w0(x,mx)=v0(x)+χ0(mx)=xMσ(q1)m+6V[∫0mx‍dth(t)2-(∫01‍dth(t)3)-1mmmmm×∫01‍dth(t)2∫0mx‍dth(t)3]=xMσ(q1)m+6Vb2m(x)b3(1)-b2(1)b3m(x)b3(1)=xMσ(q1)+Vχ0(mx).
That is a very good approximation of our exact solution (50).
It is important to notice that the choice of constantsC0,C1 was determined from the exterior boundary condition. So, we should expect the same in two-dimensional case. However the treatment of boundary conditions in two-dimensional case is much more complicated and the boundary layer is to be expected.The derived asymptotic expansion should be justified by proving the convergence. And we need the strong convergence (with corrector, of course) forwm in order to get the convergence for pm. The form of the approximation
(77)wm≈v0(x)+Vχ0(mx)+⋯.
suggests that the boundary layer phenomenon should appear on the exterior boundary ∂𝒪 since χ0 term cannot satisfy the Dirichlet condition on ∂𝒪. To get the error estimate and the strong convergence we need to handle that boundary layer. Thus, at this point we simplify the domain and the boundary condition, in order to be able to avoid it. We assume that
(78)𝒪=]0,1[×R,(79)pε(0,x2)=0,pε(1,x2)=q(x2),(80)x2↦pε(x1,x2)is1-periodic,(81)h(y1,y2)=h(y1),(82)qis1-periodic.
Now
(83)wε=Mσ(pε),withMσ(p)=∫σp‍dsμ(s)⇒wε(0,x2)=0,wε(1,x2)=Mσ(q(x2))
and x2↦wε(x1,x2) is 1-periodic.In that case we can computeχ0 and χk,k=1,2, explicitly and we can impose exterior condition on χ0. Indeed χ0 is exactly the same as in the monodimensional case; that is, it is given by (73) and (75). Obviously χ2=0 so that
(84)a22=∫01‍h(s)3ds=〈h3〉,a12=a21=0.
As for the last term
(85)χ1=-y1+(∫01‍dsh(s)3)-1∫0y1‍dsh(s)3,a11=(∫01‍dsh(s)3)-1=1〈1/h3〉.
Finally, the function v0 satisfies the boundary value problem
(86)1〈1/h3〉∂2v0∂x12+〈h3〉∂2v0∂x22=0in𝒪,v0(0,x2)=0,v0(1,x2)=Mσ(q(x2)),v0is1-periodicinx2.
It can be solved using the Fourier method, and we get
(87)v0(x1,x2)=∑k=1∞‍sh(〈h3〉〈h-3〉kπx1)×(αksinkπx2+βkcoskπx2),αk=2∫01‍Mσ(q(t))sinkπtdt,βk=2∫01‍Mσ(q(t))coskπtdt.
Since the approximation
(88)wm≈v0(x)+Vχ0(mx1)+1mχ1(mx1)∂v0∂x1(x)
now satisfies the boundary conditions on ∂𝒪, it is easy to see that
(89)|wm-(v0(x)+Vχ0(mx))|L∞(𝒪)≤C1m
follows from the maximum principle. Assuming that, for m large enough,
(90)v0(x)+Vχ0(mx)≤M+
we have
(91)|pm-H(v0(x)+Vχ0(mx))|L∞(𝒪)⟶0.
Finally
(92)pm⇀p0=∫01‍H(v0(x)+Vχ0(y1))dy1.
We have proved that.Theorem 7.
Let
p
m  be the solution to the problem (67), (79), and (80) and letv0, χ0  be defined by (87) and (42), respectively. If (90) holds, then(93)|pm-H(v0(x)+Vχ0(mx))|L∞(𝒪)⟶0,pm⇀p0=∫01‍H(v0(x)+Vχ0(y1))dy1.Remark 8.
It is important to notice thatπ0=H(p0)   satisfies
(94)∑i,k=12‍∂∂xi(aikμ(π0)∂π0∂xk)=0in𝒪,π0=qon∂𝒪,
and thus we would expect it to be the limit of pm  in analogy with the linear case. However p0≠π0.
IfV=limm→∞m|Vm|  is small, we can expand H(v0(x)+Vχ0(mx))   in powers of V and we get
(95)pm(x)≈H(v0(x)+Vχ0(mx))mmm=H(v0(x))+VH′(v0(x))χ0(y)+O(V2)mmm=π0(x)+Vμ(π0(x))χ0(y)+O(V2).
Thus
(96)p0=π0+μ(π0)〈χ0〉V+O(V2).
It can be, formally, written as
(97)pm(x)≈π0(x)+|Vm|mμ(π0(x))χ0(mx)+O(|Vm|2m2).


## 4.1. One-Dimensional Case

Ifn=1, the above problem is posed on an interval 𝒪=]0,1[. With an appropriate boundary condition
(39)pm(0)=q0,pm(1)=q1.
It forms a boundary value problem for nonlinear ODE:
(40)(hm3μ(pm)pm′)′=6Vmhm′in𝒪.
To study the asymptotic behavior of the solution with respect to m we linearize the problem using the transformation wm=Mσ(pm). To simplify, in this section we choose σ=q0 and, dropping the index σ in Mσ and Mσ+, we denote
(41)M(p)=∫q0p‍dtμ(t),M+=∫q0+∞‍dtμ(t).Theorem 5.
Let
(42)
χ
0
(
y
)
=
6
(
∫
0
1
‍
d
s
h
(
s
)
3
)
-
1
m
×
[
〈
1
h
3
〉

∫
0
y
‍
d
v
h
(
v
)
2
-

〈
1
h
2
〉

∫
0
y
‍
d
v
h
(
v
)
3
]
,
(43)
〈
·
〉
=
∫
0
1
‍
·
d
y
,
and let
H
=
M
-
1. Suppose that there exists a limit(44)V=limm→∞Vmmand that, form  large enough andM+  defined in (27), the following condition holds:(45)b3m(x)b3m(1)M(q1)+6Vmb3m(1)b2m(x)-b3m(x)b2m(1)b3m(1)≤M+,with(46)bαm(x)=∫0x‍dthm(t)α,α=2,3.Then(47)pm⇀p0=∫01‍H(xM(q1)+Vχ0(y))dymmmmmmmweak*inL∞(0,1).Proof.
Equation (40), with boundary conditions (39), can be solved by reduction to quadratures, after the substitution
(48)wm=M(pm),
with Mσ strictly increasing function defined by (12). The problem for wm now reads
(49)(hm3wm′)′=6Vmhm′in]0,1[,wm(0)=M(q0)=0,wm(1)=M(q1).
It is easy to see that (49) has a unique solution given by (25). Since M(q0)=0, (25) now reduces to
(50)wm=b3m(x)b3m(1)Mσ(q1)+6Vmb3m(1)b2m(x)-b3m(x)b2m(1)b3m(1),
where, for α=2,3,
(51)bαm(x)=∫0x‍dthm(t)α=1m∫0mx‍dsh(s)α⟶x∫01‍dsh(s)α,mmmmmmmmmasm⟶+∞.
Now pm, the solution to the problem (40), exists if (45) is fulfilled. The second term in (50) thus obviously tends to xM(q1), as m→+∞. The last term is more interesting. The denominator tends to
(52)〈1h3〉=∫01‍dsh(s)3.
As for its numerator, we have
(53)Vm[b3m(1)b2m(x)-b3m(x)b2m(1)]mm=Vm[∫01‍∫0x‍(1hm(t)3hm(s)2-1hm(t)2h(s)m3)dsdt]mm=Vmm(〈1h3〉∫0mx‍dvh(v)2-〈1h2〉∫0mx‍dvh(v)3).
Suppose that(54)V=limm→+∞Vmm
and denote
(55)G(y)=〈1h3〉∫0y‍dvh(v)2-〈1h2〉∫0y‍dvh(v)3.
Obviously the function G is periodic with period 1 so that, due to the standard periodicity lemma (see, e.g., [19]), as m→+∞,
(56)G(mx)⇀〈G〉=∫01‍G(y)dy,weak*inL∞(0,1).
By direct computation(57)〈G〉=〈yh3〉〈1h2〉-〈1h3〉〈yh2〉=∫∫01s-th(s)3h(t)2dsdt.
Thus
(58)Vm[b3m(1)b2m(x)-b3m(x)b2m(1)]⇀V〈G〉.
Now, denoting(59)χ0(y)=G(y)b3(1)
we have
(60)wm≈xM(q1)+Vχ0(mx),
and thus
(61)wm⇀w0=xM(q1)+V〈χ0〉weak*inL∞(0,1).
However, we are not interested in convergence of the auxiliary functionwm but in the convergence of the pressure pm. Since (45) is assumed to be true, we can define pm=H(wm), where H=M-1, and we have
(62)pm⇀p0=∫01‍H(xM(q1)+Vχ0(y))dymmmmmweak*inL∞(0,1).Deriving the expression on the right-hand side, we obtain the effective pressure drop in the form(63)p0′(x)=M(q1)∫01‍μ(xM(q1)+Vχ0(y))dy(64)=∫q0q1‍dτμ(τ)∫01‍μ(xM(q1)+Vχ0(y))dy=∫q0q1‍dτμ(τ)∫01‍μ(x∫q0q1‍dsμ(s)+Vχ0(y))dy.
As we can see, the pressure drop is not constant, as for the Newtonian flow. The interesting effect appears if V≠0 because, in that case, the expressions for the pressure and for the pressure drop are nonlocal due to the integral with respect to y. That phenomenon is entirely due to the fact that the viscosity is depending on the pressure.


## 4.2. Two-Dimensional Case

We suppose here that the functionhm is constructed from positive, smooth, Y-periodic function h:R2→[d0,+∞[, d0>0, Y=]0,1[2 in the same way as before, that is, by taking
(65)hm(x)=h(mx).
We have seen in the previous section that interesting effects happen only if we assume that
(66)Vm=mV.
In that case our equation reads
(67)div(hm3μ(pm)∇pm)=6Vm·∇hmin𝒪.
As we did in the existence analysis and in the previous section, we linearize the equation by substitution
(68)wm=Mσ(pm),
where the function Mσ is defined by (12). Now wm satisfies
(69)div(h3M∇wm)=6Vm·∇hmin𝒪.We postulate the asymptotic expansion in the form(70)wm(x)≈w0(x,y)+1mw1(x,y)m+1m2w2(x,y)+⋯,y=mx.
All functions are assumed to be Y-periodic in y variable.Plugging that in (69) and collecting, formally, terms with equal powers of m, we getm2: divy(h3∇yw0)=6V·∇yhm: divy(h3∇yw1)+divy(h3∇xw0)+h3divx∇yw0=01:divy(h3∇yw2)+divy(h3∇xw1)+h3divx∇yw1+h3Δxw0=0.Denoting(71)V=|V|,V0=VV
we have
(72)w0(x,y)=v0(x)+Vχ0(y),divy(h3∇yχ0)=6V0·∇yhinY,w1(x,y)=∑k=12‍χk(y)∂v0∂xk(x)+v1(x),divy(h3∇y(χk+yk))=0,∑k,ℓ=12‍akℓ∂2v0∂xk∂xℓ=0in𝒪,akℓ=∫Y‍h3∂∂yℓ(χk+yk)dy.Remark 6.
The same computation can be done in one-dimensional case and it gives(73)(h3χ0′)′=6h′⇒χ0=6∫0y‍dth(t)2+C0∫0y‍dth(t)3+C1,(74)v0=xM(q1).
Constants C0,C1 are chosen in a way that boundary conditions χ0(0)=χ0(1)=0  are met, and it follows that
(75)C1=0,C0=-6(∫01‍dth(t)3)-1∫01‍dth(t)2.
Then
(76)wm≈w0(x,mx)=v0(x)+χ0(mx)=xMσ(q1)m+6V[∫0mx‍dth(t)2-(∫01‍dth(t)3)-1mmmmm×∫01‍dth(t)2∫0mx‍dth(t)3]=xMσ(q1)m+6Vb2m(x)b3(1)-b2(1)b3m(x)b3(1)=xMσ(q1)+Vχ0(mx).
That is a very good approximation of our exact solution (50).
It is important to notice that the choice of constantsC0,C1 was determined from the exterior boundary condition. So, we should expect the same in two-dimensional case. However the treatment of boundary conditions in two-dimensional case is much more complicated and the boundary layer is to be expected.The derived asymptotic expansion should be justified by proving the convergence. And we need the strong convergence (with corrector, of course) forwm in order to get the convergence for pm. The form of the approximation
(77)wm≈v0(x)+Vχ0(mx)+⋯.
suggests that the boundary layer phenomenon should appear on the exterior boundary ∂𝒪 since χ0 term cannot satisfy the Dirichlet condition on ∂𝒪. To get the error estimate and the strong convergence we need to handle that boundary layer. Thus, at this point we simplify the domain and the boundary condition, in order to be able to avoid it. We assume that
(78)𝒪=]0,1[×R,(79)pε(0,x2)=0,pε(1,x2)=q(x2),(80)x2↦pε(x1,x2)is1-periodic,(81)h(y1,y2)=h(y1),(82)qis1-periodic.
Now
(83)wε=Mσ(pε),withMσ(p)=∫σp‍dsμ(s)⇒wε(0,x2)=0,wε(1,x2)=Mσ(q(x2))
and x2↦wε(x1,x2) is 1-periodic.In that case we can computeχ0 and χk,k=1,2, explicitly and we can impose exterior condition on χ0. Indeed χ0 is exactly the same as in the monodimensional case; that is, it is given by (73) and (75). Obviously χ2=0 so that
(84)a22=∫01‍h(s)3ds=〈h3〉,a12=a21=0.
As for the last term
(85)χ1=-y1+(∫01‍dsh(s)3)-1∫0y1‍dsh(s)3,a11=(∫01‍dsh(s)3)-1=1〈1/h3〉.
Finally, the function v0 satisfies the boundary value problem
(86)1〈1/h3〉∂2v0∂x12+〈h3〉∂2v0∂x22=0in𝒪,v0(0,x2)=0,v0(1,x2)=Mσ(q(x2)),v0is1-periodicinx2.
It can be solved using the Fourier method, and we get
(87)v0(x1,x2)=∑k=1∞‍sh(〈h3〉〈h-3〉kπx1)×(αksinkπx2+βkcoskπx2),αk=2∫01‍Mσ(q(t))sinkπtdt,βk=2∫01‍Mσ(q(t))coskπtdt.
Since the approximation
(88)wm≈v0(x)+Vχ0(mx1)+1mχ1(mx1)∂v0∂x1(x)
now satisfies the boundary conditions on ∂𝒪, it is easy to see that
(89)|wm-(v0(x)+Vχ0(mx))|L∞(𝒪)≤C1m
follows from the maximum principle. Assuming that, for m large enough,
(90)v0(x)+Vχ0(mx)≤M+
we have
(91)|pm-H(v0(x)+Vχ0(mx))|L∞(𝒪)⟶0.
Finally
(92)pm⇀p0=∫01‍H(v0(x)+Vχ0(y1))dy1.
We have proved that.Theorem 7.
Let
p
m  be the solution to the problem (67), (79), and (80) and letv0, χ0  be defined by (87) and (42), respectively. If (90) holds, then(93)|pm-H(v0(x)+Vχ0(mx))|L∞(𝒪)⟶0,pm⇀p0=∫01‍H(v0(x)+Vχ0(y1))dy1.Remark 8.
It is important to notice thatπ0=H(p0)   satisfies
(94)∑i,k=12‍∂∂xi(aikμ(π0)∂π0∂xk)=0in𝒪,π0=qon∂𝒪,
and thus we would expect it to be the limit of pm  in analogy with the linear case. However p0≠π0.
IfV=limm→∞m|Vm|  is small, we can expand H(v0(x)+Vχ0(mx))   in powers of V and we get
(95)pm(x)≈H(v0(x)+Vχ0(mx))mmm=H(v0(x))+VH′(v0(x))χ0(y)+O(V2)mmm=π0(x)+Vμ(π0(x))χ0(y)+O(V2).
Thus
(96)p0=π0+μ(π0)〈χ0〉V+O(V2).
It can be, formally, written as
(97)pm(x)≈π0(x)+|Vm|mμ(π0(x))χ0(mx)+O(|Vm|2m2).



---

*Source: 710214-2013-04-11.xml*