## Abstract

Lubrication theory is broadly applicable to the flow characterization of thin fluid films and the motion of particles near surfaces. We offer an extension to lubrication theory by starting with Stokes equations and considering higher-order terms in a systematic perturbation expansion to describe the fluid flow in a channel with features of a modest aspect ratio. Experimental results qualitatively confirm the higher-order analytical solutions, while numerical results are in very good agreement with the higher-order analytical results. We show that the extended lubrication theory is a robust tool for an accurate estimate of pressure drop in channels with shape changes on the order of the channel height, accounting for both smooth and sharp changes in geometry.

## 1. Introduction

Lubrication theory is an approximation to the Navier–Stokes and continuity equations at low Reynolds numbers for narrow geometries with slow changes in curvature [1–3]. The approach is used regularly to describe the velocity field and pressure gradient in fluid film lubricants [4,5], the motion of particles within a fluid and near boundaries [6,7], the fluid flow passing through a microchannel with a known geometry [8–11], flow driven by the contracting walls of a soft channel [12], e.g. an insect’s trachea [13,14], and the flow of thin liquid films with free surfaces [15,16], e.g. when a droplet wets a solid surface [17,18] or in spin coating over topographically patterned surfaces [19,20]. Classical lubrication theory (CLT) is suitable for all of the above cases provided that the typical magnitude of the boundary slope is sufficiently small, typically on the order of

In this study, we obtain higher-order terms of the lubrication approximation and present an extension to lubrication theory, which we refer to as extended lubrication theory (ELT), to address two limitations of CLT. First, the use of ELT is no longer limited to small gaps and thin films. Second, the boundaries can be described by any mathematical shape function with modest curvatures as long as they are continuous and differentiable. In addition, we show how the differentiability condition may be relaxed at low Reynolds numbers, at least in practice, by considering geometries that are piecewise differentiable. We compare the results of different orders of the analytical solutions with experimental results and also with direct numerical solutions of the Navier–Stokes equations to define a threshold for considering higher-order terms in the solution. The theory will break down for sufficiently large wall slopes and then if an analysis is desired, an approach with matched asymptotic expansions is necessary, e.g. Cormack *et al.* [23]. Note that a related idea emphasizing the first correction to the classical lubrication approach, but focused on the mathematical form of the differential equations, is given by Marušić-Paloka *et al.* [24].

Our approach here can be used to describe the flow in a channel with a large constriction made by a trapped cell, particle, bubble or droplet, or as a result of buckled, crumpled or swollen walls. We validate this approach for large non-uniformities within the fluid channel by varying two characteristic geometric parameters, the constriction’s length and amplitude, from zero to the channel height, e.g. the corresponding dimensionless geometric parameters vary from zero to one. We note that in addition to numerical solutions, two types of analytical calculations have been considered in the literature. One method is domain perturbation [25–29] where the constriction amplitude is much smaller than the channel height. The second method involves corrections to leading-order lubrication theory [10,30–32] while the constriction amplitude can be as high as the channel height. We show the applicability of our approach to different geometries and provide the accuracy of such results as two non-dimensional geometric parameters along with the Reynolds number vary. Finally, we note that a similar higher-order expansion was used recently in a study of an electrokinetic flow in a channel of non-uniform shape [33].

## 2. Theoretical approximation

We consider working to higher order in traditional lubrication theory to describe fluid flow in non-uniform channel shapes with modest aspect ratios. Thus, we consider incompressible, steady, two-dimensional pressure-driven flow in a channel with shape *y*=*h*(*x*)=*h*_{0}*H*(*X*), where *X*=*x*/*L*_{0}, *L*_{0} is the channel length, *h*_{0} is a characteristic channel height, *H*(*X*) is a normalized shape function and *δ*=*h*_{0}/*L*_{0}≪1. A typical geometry in the form of a constriction is shown in figure 1*a*. We assume that the Reynolds number is small and so consider the continuity and Stokes equations
**u**=(*u*,*v*) is the velocity field and *μ* is the fluid viscosity. We denote the constant flow rate (per unit width) as *q*_{0}. A discussion of the flow based on the full Navier–Stokes equations is given in §4.

Consistent with the traditional lubrication approximation, we analyse equation (2.1) and choose to introduce dimensionless variables according to
*v* is *u*, which is the component of flow along the channel axis. The dimensionless equations corresponding to (2.1) for a two-dimensional flow are
*a*
*b*
*c*These equations are to be solved with boundary conditions
*δ*^{2}. However, we later introduce another dimensionless parameter, λ, coming from the shape function and corresponding to the normalized constriction amplitude. Also, we will show below that, for sinusoidal shape variations with dimensionless amplitude λ (figure 1*b*), the effective small parameter is (*δ*λ)^{2}, corresponding to the square of the typical wall slope.

### (a) Perturbation expansion and the leading-order result from lubrication theory

Our first steps follow standard discussions in textbooks, e.g. [2]. Because the problem statement only involves *δ*^{2}, which is assumed to be small, we seek a solution to (2.3) of the form
*a*
*b*
*c*

At leading order, we have the familiar classical lubrication problem
*a*
*b*
*c*with *U*_{0}=0 at *Y* =0 and *H*(*X*). The solution is
*X*, follows from applying the integral constraint

To provide an example, we consider the shape function
*b*. Here, λ corresponds to the dimensionless amplitude of the change in the channel height. Thus, the wall slope has magnitude |d*h*/d*x*|=*O*(*h*_{0}/*L*_{0}λ)=*O*(*δ*λ).

The leading-order pressure drop Δ*P*_{0} is then calculated to be (the integration was accomplished with Mathematica)

Before proceeding further, we determine the velocity component *V*_{0}(*X*,*Y*) using the continuity equation. Although equation (2.6a) is first order in *Y* , we expect it to satisfy two boundary conditions, as *V*_{0}(*X*,0)=*V*_{0}(*X*,*H*(*X*))=0. Using the continuity equation, and imposing *V*_{0}(*X*,0)=0, we have
*X* derivatives. We then note that, at *Y* =*H*(*X*), direct differentiation shows that (2.13) yields *V*_{0}(*X*,*H*(*X*))=0. Alternatively, we can write
*Y* -component of velocity at every order in the analysis below and the no-slip boundary condition is satisfied for both velocity components of **u**.

### (b) The O ( δ 2 ) term in the perturbation expansion

In most calculations using lubrication theory, the development is truncated with the leading-order term calculated in the preceding section, which is well-known from countless applications. Here, our interest is to improve the approximation by including additional terms in the perturbation solution. At the next order, *a*
*b*
*c*with boundary conditions *U*_{2}=0 at *Y* =0 and *H*(*X*), and *P*_{2}=*P*_{2}(−1)−*P*_{2}(1) needed to enforce the constraint on the flux.

We can integrate the last equation of (2.15), and use continuity, to obtain
*c*_{3}(*X*) is allowed by the integration. With this pressure distribution, we use the *X*-momentum equation to find
*U*_{0} is given by equation (2.9). Upon integration, and application of the boundary conditions, we find
*X*.

As we have accounted for the specified dimensionless flow rate at leading order, we now require *X*-component of the velocity *U*_{2}(*X*,*Y*) for any shape function *H*(*X*). To continue with the example of figure 1, we again use the shape function in equation (2.10). Integrating (2.16), taking into account that ∂*U*_{0}/∂*X* vanishes as *P*_{2} at this order as:

We determine *V*_{2}(*X*,*Y*) using the continuity equation and imposing *V*_{2}(*X*,0)=0, which leads to the expression

This equation only involves the shape function *H*(*X*) because d*c*_{3}/d*X* is given in (2.19). As in the previous section, it can be verified that *V*_{2}(*X*,*H*(*X*))=0.

We conclude this subsection by noting that we have now established the net pressure drop at this order:
*O*((*δ*λ)^{2}), which corresponds to the square of the wall slope.

### (c) The perturbation expansion at O ( δ 4 )

It is useful to go one step further simply to illustrate that the basic analytical steps carry through at every order. The higher-order terms help to provide a better representation of flows in geometries with more rapid shape variations. We can continue these basic steps at *a*
*b*
*c*with *U*_{4}=0 at *Y* =0 and *H*(*X*), and *Y* -momentum equation can be integrated, which, after using the continuity equation, yields
*X*-momentum equation we have
*U*_{2}(*X*,*Y*) is known from equation (2.18), then we calculate
*U*_{4}=0 at *Y* =0 and *H*(*X*) to arrive at
*c*_{5}/d*X*:

Equations (2.19), (2.28) and (2.29) give the *X*-component velocity at this order for any choice of the shape function *H*(*X*). We determine the correction to the pressure drop as:

For a given flow rate (*q*_{0}), we have determined the dimensionless pressure drop *δ* and λ, where Δ*p*_{ measured} is the difference in pressure measured at the two ends of the constriction. In particular,

We next describe experiments and numerical simulations to confirm the improved description offered by these additional terms in the lubrication approximation.

## 3. Experimental verification

Our experimental set-up consists of a long channel (200 mm) with a rectangular cross section (5 ×5 mm) and an obstruction in the middle (figure 2*a*), with *L*_{0}≈*h*_{0}. The sides of the channel were cut from acrylic sheets (8560K211, McMaster-Carr) using a laser cutter (Epilog Mini Laser 24, 60 W) and bonded together using an acrylic capillary cement (10705, TAP Plastics). We varied the arch size and shape, similar to figure 1, while keeping all other geometrical parameters constant between different tests. We recognize that our theory is two dimensional and the experimental geometry is three dimensional. However, as the arch amplitude increases, the flow through the narrow gap better approximates a two-dimensional flow, which allows us to highlight the impact on the pressure drop of systematically increasing the magnitude of the shape variations.

The pressure drop within the channel between two fixed points, symmetrically located on each side of the arch, was measured using a sensitive differential pressure sensor (CPCL04D, Honeywell). By keeping the flow rate constant for different tests, the pressure drop across the arch was then obtained by subtracting the pressure drop within the flat part of the channel from the total pressure drop between those fixed points. The fluid was chosen to be a standard viscosity oil (N1000, Cannon Instrument) and the temperature was kept at 22.5±0.5^{°}C, resulting in a viscosity of 2.45 Pa s and density of 848 kg m^{−3}. We used a syringe pump (PHD Ultra CP, Harvard Apparatus) to apply a fixed flow at a rate of 14.4 mm^{3} s^{−1} so that

We repeated each test at least three times and the average measured pressure drops along with their associated standard deviations are shown in figure 2*b*, where *δ*≈1 and λ is varied from zero (flow in a straight channel) to λ≈0.9. We note that higher orders of the analytical solutions are in better agreement with the experimental results, while the CLT (red dotted line) underestimates the results by about 40%. We also note that the analytical solutions are obtained in a two-dimensional channel, while the experimental results are from a three-dimensional channel, and this is one of the reasons for the difference between the analytical and experimental results. Although the experimental results confirm the trend in the analytical solutions, it may not be feasible to conduct experiments for every case to verify the higher orders of the analytical solutions due to the sensitivity and the difficulty of such experiments. So we performed numerical simulations for a variety of cases and compare those results with the analytical solutions in the next section.

## 4. Numerical simulations

As an additional route to test the analytical results above, we seek to numerically solve the Navier–Stokes equation in their full form for incompressible, steady two-dimensional flow. Again, we assume a specified flow rate and seek the corresponding pressure drop as geometric parameters are varied. The same scalings and dimensionless parameters introduced in equation (2.2) are used to obtain the dimensionless continuity and Navier–Stokes equations
*a*
*b*
*c*where *ρ* is the fluid density and *μ* is the fluid viscosity. In the lubrication literature, it is common to define the reduced Reynolds number, *δ* and

A numerical solver often uses the weak form of (4.1). To do so, we consider an arbitrary pair of *P* and **U**=(*U*,*V*) to be a solution to the dimensionless continuity and Navier–Stokes equations (4.1) for a steady and incompressible flow. If these equations are multiplied by any pressure and velocity basis functions, i.e. (*q*,*ν*_{1},*ν*_{2}), and integrated over the domain *Ω*, the pair is still a solution, and satisfies the new equations. We then reduce all the second-order terms to first-order ones using Gauss’ theorem and neglect the boundary integrals as they are usually handled separately in finite element packages. Therefore, we have the following weak form of the continuity and Navier–Stokes equations that can directly be used within numerical solvers:
*a*
*b*
*c*

We used COMSOL 4.3b mainly as a PDE solver to the above equations. These dimensionless equations depend on two dimensionless parameters: Reynolds number and the geometric variable *δ*=*h*_{0}/*L*_{0}. In addition, the shape function (2.10) adds another dimensionless parameter, λ, to this system. So we performed the numerical simulations first for *δ* from 0.2 to 1. The dimensionless channel height was *H*_{0}=1 and we applied a flow at the inlet at a fixed rate with a parabolic velocity profile of 6(*Y* −*Y* ^{2}) so that *Q*=1. The outlet pressure was also set to zero. After solving the governing equations, the pressure drop was measured at two cross sections that are symmetrically located on each side of the arch and separated by 2*L*_{0}. These scales respect the non-dimensionalization in §2. While we recognize that our theory, which is based on Stokes equations, is strictly valid only when *h*_{0} and *L*_{0}), the ratio of the inertia to the viscous terms in the Navier–Stokes equations involves the product of *δ*. Thus, we also performed some numerical calculations with finite

The comparison between simulation results and different orders of the analytical solutions for a channel with a shape function provided in (2.10) is shown in figure 4. As expected from textbook discussions of CLT (see also [29]), when *δ* is small, i.e. *δ*≤0.2, the CLT as well as higher-order ELTs estimate the pressure drop accurately (the errors are within 5% when *δ* (including term *a*). We also note that if a maximum 4% error is acceptable, the second-order ELT would be adequate to estimate the pressure drop for a simple shape like (2.10) while *δ*≤1 and *b*), even the higher orders of the analytical solutions do not estimate the pressure drop accurately as inertial forces become more dominant than the viscous forces and can no longer be ignored. As figure 4*b* shows this might be considered as an upper limit to rely on such analytical solutions.

We then investigate the use of higher-order ELT in applications where a channel, instead of an obstruction, has a convex shape. We choose the shape function (2.10) while varying λ from 0 (no bulge) to −1 (bulge with the size of the channel height) (figure 5*a*). This convexity alters the flow profile (figure 5*b*) and reduces the pressure drop within the channel. We performed numerical analyses for different *c*. As we have shown that even for *δ*=1 the ELT provides a reasonable approximation to the full numerical simulations (figure 4*a*), here we choose *δ*=1. For a channel with a sharp convexity, only fourth-order ELT and higher may accurately estimate the pressure drop (within 5% error), while CLT and the second-order ELT estimations differ from the simulation results by about 30% and 20%, respectively. Therefore, higher orders of the analytical solutions significantly improve the estimation of pressure drop within a channel with a significant change in geometry.

Until now, we have used shape functions that are entirely differentiable. This condition may not be met in models of all applications. Here, we provide of an example showing that ELT can be applied to a channel whose shape is continuous but piecewise differentiable, i.e. the shape function may have finite non-differentiable points. A channel with such a shape function can be divided into smaller sections where each part has a differentiable shape. The pressure drop within each piece is estimated using ELT. Since the flow is laminar, the total pressure drop across the original channel is the summation of the pressure drops within each part. For example, consider a shape function with a single non-differentiable point as follows:
*γ*<1 is a dimensionless parameter that determines the location of the discontinuity in slope and 1−λ gives the minimum gap height. This shape function is plotted for *γ*=0.75 in figure 6*a*. Following the same procedure introduced in §2, the pressure drop is
*P*_{0}=12 (2−λ)/(1−λ)^{2}, and the terms inside the parentheses correspond to CLT, second-order ELT, fourth-order ELT and so on, respectively. We used the same shape function (4.3) to perform numerical simulations, and the comparison is shown in figure 6*b*. For channels with small *δ* (*δ*≤0.2), theoretical and numerical results are in good agreement, while for channels with *δ*≈1, the pressure drop estimated using the higher orders of the analytical solutions follow the numerical results more closely. These results show that ELT can be applied to a channel with a piecewise differentiable shape function.

We note that shape functions can appear on both sides of a channel and with appropriate boundary conditions (2.4), the same procedure can be followed to find analytical solutions at different orders. This method can potentially be extended to derive analytical solutions for channels with two-dimensional shape functions on one or both sides. We note that, for some extensions, the algebra becomes more cumbersome and intermediate integrals may need to be accomplished numerically, which may partly defeat the utility of the current approach.

## 5. Conclusion

We extended the lubrication approximation by obtaining higher-order terms in a systematic perturbation analysis and compared the analytical results with experiment and numerical simulations. Experimental results were closer to higher-order analytical solutions when the gap was narrow so that the two-dimensional approximation was appropriate. Very good agreement was found between higher-order analytical solutions and the simulation results, confirming that, for channels with a high aspect ratio, the higher-order terms of the ELT results in a significant improvement in accuracy when compared with the CLT. Moreover, though domain perturbations are unsuccessful when the shape changes are modest, ELT has been shown to be quite successful, at least for the channel shapes studied here. Nevertheless, if the wall slope is too steep, we can expect this lubrication approach to fail. Finally, for low Reynolds numbers, simple piecewise differentiable shape functions can be used with the analytical solutions obtained in this study, which provides a robust tool to accurately estimate the pressure drop in a channel with positive or negative constrictions, whose changes in height are comparable with its length.

## Data accessibility

All datasets presented in this work would be available upon request.

## Authors' contributions

B.T., G.F., D.P.H. and H.A.S. derived the analytical approximations. B.T. and D.P.H. carried out the experiments. B.T. performed the numerical analyses. B.T., G.F., D.P.H. and H.A.S. wrote the paper.

## Competing interests

The authors have no competing interests.

## Funding

This work was supported by the National Science Foundation (B.T. and D.P.H., CMMI-1505125).

- Received March 31, 2017.
- Accepted September 5, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.