2 Lubrication approximation for shallow deformable channels

We consider a channel of length $\ell$ , width $w$ and height $h$ , where $h\ll w\ll \ell$ . The upper wall of the channel is soft and deformable, as is the case when a rigid channel is sealed by a thin elastic film. The flow is in the positive $z$ -direction. Due to the normal stresses from the flow on the walls, the (soft) top wall of the channel deforms out of the $(x,z)$ -plane in the positive $y$ -direction, so that the steady shape of the channel’s top wall is given by $y=h(x,z)=h_{0}+u(x,z)$ , as shown in figure 1. For the moment, we make no assumption on the magnitude of the displacement, however, we expect that $u>0$ and $u\ll w$ for the types of problems of interest herein.

Figure 1.

Three-dimensional schematic of the deformed geometry of a channel of length $\ell$ with an initially rectangular cross-section ( $w\times h_{0}$ ). The volumetric flow rate is denoted by $q$ and $u(x,z)$ is the top wall’s deformation (the remaining walls are assumed rigid). Although not shown in this schematic, the top wall has a thickness $t$ .

Now, consider the incompressible flow of a Newtonian fluid of density $\unicode[STIX]{x1D71A}$ and viscosity $\unicode[STIX]{x1D707}$ . As is typical for microfluidic devices (Stone, Stroock & Ajdari Reference Stone, Stroock and Ajdari2004; Squires & Quake Reference Squires and Quake2005; Bruus Reference Bruus2008), we are interested in the limit in which the effective Reynolds number is small, i.e. $\unicode[STIX]{x1D716}Re=(h_{0}/\ell )\unicode[STIX]{x1D71A}{\mathcal{V}}_{c}h_{0}/\unicode[STIX]{x1D707}\ll 1$ , where ${\mathcal{V}}_{c}$ is a characteristic axial velocity and $\unicode[STIX]{x1D716}=h_{0}/\ell \ll 1$ . In this limit, fluid inertia is negligible compared to viscous stresses, and the flow is governed by the Stokes equations (Happel & Brenner Reference Happel and Brenner1983; Bruus Reference Bruus2008):

(2.1a,b ) $$\begin{eqnarray}\mathbf{0}=-\unicode[STIX]{x1D735}p+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{v},\quad \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0,\end{eqnarray}$$

where $\boldsymbol{v}=(v_{x},v_{y},v_{z})$ and body forces have been neglected. The flow is also subject to the no-slip and no-penetration boundary conditions along the walls of the channel:

(2.2) $$\begin{eqnarray}\boldsymbol{v}=\mathbf{0}\quad \text{at }\left\{\begin{array}{@{}l@{}}y=0,h(x,z),\quad \\ x=\pm w/2.\quad \end{array}\right.\end{eqnarray}$$

Note that the velocity vanishes along $y=h$ , rather than being $v_{x}=\text{d}h/\text{d}t$ , etc. because we have already assumed the fluid flow and the structural deformation to be steady (i.e. independent of time).

Let us introduce the dimensionless variables

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}x=wX,\quad y=h_{0}Y,\quad z=\ell Z,\quad h(x,z)=h_{0}H(X,Z),\quad u(x,z)=u_{\mathit{c}}U(X,Z),\\ v_{x}(x,y,z)={\displaystyle \frac{\unicode[STIX]{x1D716}}{\unicode[STIX]{x1D6FF}}}{\mathcal{V}}_{\mathit{c}}V_{X}(X,Y,Z),\quad v_{y}(x,y,z)=\unicode[STIX]{x1D716}{\mathcal{V}}_{\mathit{c}}V_{Y}(X,Y,Z),\\ v_{z}(x,y,z)={\mathcal{V}}_{\mathit{c}}V_{Z}(X,Y,Z),\quad p(x,y,z)=\unicode[STIX]{x0394}p\,P(X,Y,Z),\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x0394}p$ , $u_{\mathit{c}}$ and ${\mathcal{V}}_{\mathit{c}}$ are a characteristic pressure scale, a characteristic top wall deformation scale and a characteristic axial velocity scale, respectively; we have also defined $\unicode[STIX]{x1D6FF}:=h_{0}/w$ and $\unicode[STIX]{x1D716}:=h_{0}/\ell$ . The characteristic top wall deformation $u_{\mathit{c}}$ is not an independent parameter and will be determined below, upon scaling the governing equations. First, consider the continuity equation (2.1b ), which becomes

(2.4) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}V_{X}}{\unicode[STIX]{x2202}X}+\frac{\unicode[STIX]{x2202}V_{Y}}{\unicode[STIX]{x2202}Y}+\frac{\unicode[STIX]{x2202}V_{Z}}{\unicode[STIX]{x2202}Z}=0.\end{eqnarray}$$

By convention, we sought a balance in the continuity equation, which is how the velocity scales in (2.3) were set.

A microchannel can be operated in one of two regimes: either the pressure drop $\unicode[STIX]{x0394}p=p(0)-p(\ell )$ is imposed, where we assume that the outlet is always open to the atmosphere so we set $p(\ell )=0$ (then, the inlet pressure is $p(0)=\unicode[STIX]{x0394}p$ ), or the flow rate $q$ is imposed, which is equivalent to imposing the mean flow velocity $\langle v_{z}\rangle \propto q/(h_{0}w)$ . In both cases, $\unicode[STIX]{x0394}p$ and ${\mathcal{V}}_{\mathit{c}}$ are not independent scales but coupled by the flow physics. In the first case, the characteristic velocity scale is set by pressure drop $\unicode[STIX]{x0394}p$ : ${\mathcal{V}}_{\mathit{c}}=h_{0}^{2}\unicode[STIX]{x0394}p/(\unicode[STIX]{x1D707}\ell )$ ; in the second case, the pressure drop is set by the characteristic axial velocity scale ${\mathcal{V}}_{\mathit{c}}=\langle v_{z}\rangle$ : $\unicode[STIX]{x0394}p=\unicode[STIX]{x1D707}\langle v_{z}\rangle \ell /h_{0}^{2}$ . In either case, equation (2.1a ) adopts the same form:

(2.5a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FF}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{X}}{\unicode[STIX]{x2202}X^{2}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{X}}{\unicode[STIX]{x2202}Y^{2}}+\unicode[STIX]{x1D716}^{4}\frac{\unicode[STIX]{x2202}^{2}V_{X}}{\unicode[STIX]{x2202}Z^{2}} & = & \displaystyle \unicode[STIX]{x1D6FF}^{2}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X},\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FF}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{Y}}{\unicode[STIX]{x2202}X^{2}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{Y}}{\unicode[STIX]{x2202}Y^{2}}+\unicode[STIX]{x1D716}^{4}\frac{\unicode[STIX]{x2202}^{2}V_{Y}}{\unicode[STIX]{x2202}Z^{2}} & = & \displaystyle \frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Y},\end{eqnarray}$$

(2.5c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FF}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{Z}}{\unicode[STIX]{x2202}X^{2}}+\frac{\unicode[STIX]{x2202}^{2}V_{Z}}{\unicode[STIX]{x2202}Y^{2}}+\unicode[STIX]{x1D716}^{2}\frac{\unicode[STIX]{x2202}^{2}V_{Z}}{\unicode[STIX]{x2202}Z^{2}} & = & \displaystyle \frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Z}.\end{eqnarray}$$

Above, we assumed that the channel is long and thin, which leads to a natural ordering of the small (dimensionless) parameters:

(2.6) $$\begin{eqnarray}0<\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D6FF}\ll 1.\end{eqnarray}$$

We are interested in the leading-order asymptotic behaviour under the ordering (2.6). From (2.5), we have

(2.7a-c ) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{2}V_{Z}}{\unicode[STIX]{x2202}Y^{2}}=\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Z},\quad 0=\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Y},\quad 0=\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}X},\end{eqnarray}$$

which is the familiar lubrication approximation. From (2.7b,c ), we first infer that the leading-order pressure does not depend on the cross-sectional coordinates $(X,Y)$ . Then, after integrating (2.7a ) and enforcing the no-slip boundary condition from (2.2), i.e. $V_{Z}(X,0,Z)=V_{Z}(X,H(X,Z),Z)=0$ $\forall X,Z$ , the ‘standard’ lubrication theory result for the axial velocity follows:

(2.8) $$\begin{eqnarray}V_{Z}(X,Y,Z)=-\frac{1}{2}\frac{\text{d}P}{\text{d}Z}(H(X,Z)-Y)Y.\end{eqnarray}$$

We remind the reader that $\text{d}P/\text{d}Z<0$ in our convention. This leading-order axial velocity does not satisfy the no-slip boundary condition at the side walls ( $X=\pm 1/2$ ). Appendix B discusses the effect of the lateral side walls, showing their effect is indeed small for $\unicode[STIX]{x1D6FF}\ll 1$ .

From the velocity scalings given in (2.3), it is clear that the cross-sectional velocity components are much smaller than the axial one. Thus, to the leading order in the two small parameters, the velocity field in this lubrication approximation is dominated by the axial component $V_{Z}$ . Terms beyond the leading-order ones can be computed in the standard way, by positing a regular perturbation expansion in the two small parameters $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6FF}$ .

Next, we wish to relate the volumetric flow rate though the channel,

(2.9) $$\begin{eqnarray}Q:=\int _{-1/2}^{+1/2}\int _{0}^{H(X,Z)}V_{Z}(X,Y,Z)\,\text{d}Y\,\text{d}X,\end{eqnarray}$$

to the streamwise gradient of the pressure $\text{d}P/\text{d}Z$ . Since the flow is steady, $Q$ is constant. Then, substituting $V_{Z}$ from (2.8) into (2.9) and integrating over $Y$ , we get (to the leading order in $\unicode[STIX]{x1D6FF}$ and $\unicode[STIX]{x1D716}$ ) the ‘standard’ lubrication theory relation

(2.10) $$\begin{eqnarray}Q=-\frac{1}{12}\frac{\text{d}P}{\text{d}Z}\int _{-1/2}^{+1/2}H(X,Z)^{3}\,\text{d}X.\end{eqnarray}$$

Thus, for a given shape of the top wall of the channel, we can find $Q$ from (2.10). Within the lubrication approximation, the first correction in $\unicode[STIX]{x1D6FF}$ to (2.10) can be computed as shown in appendix B.

3 Shape of the deformed channel

In terms of the dimensionless variables from (2.3), the channel shape can thus be expressed as

(3.1) $$\begin{eqnarray}H(X,Z)=1+\unicode[STIX]{x1D6FD}U(X,Z),\end{eqnarray}$$

where we have set $\unicode[STIX]{x1D6FD}:=u_{\mathit{c}}/h_{0}$ . Here, the dimensionless group $\unicode[STIX]{x1D6FD}$ controls the compliance of the top wall: for $\unicode[STIX]{x1D6FD}\ll 1$ , the top wall is stiff (equivalently, its deformation is small compared to the undeformed height), and for $\unicode[STIX]{x1D6FD}\gg 1$ it is soft (equivalently, its deformation is large compared to the undeformed height). Within the lubrication approximation, we do not need to make any assumptions on the smallness of $\unicode[STIX]{x1D6FD}$ at this stage, however, we expect $\unicode[STIX]{x1D6FD}>0$ and $\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D6FF}\ll 1$ as stated above in § 2 (see also appendix B).

If the maximum displacement $\max _{x,z}|u(x,z)|$ of the top wall can be assumed small compared to its thickness $t$ , and its thickness $t$ is small compared to its width $w$ , then, from the theory of linear elasticity, we know that the steady-state displacement $u(x,z)$ satisfies the Kirchhoff–Love equation for isotropic quasi-static bending of a plate under a transverse load due to the fluid pressure (Love Reference Love1888; Timoshenko & Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1959; Landau & Lifshitz Reference Landau and Lifshitz1986):

(3.2) $$\begin{eqnarray}B\unicode[STIX]{x1D6FB}_{\Vert }^{2}\unicode[STIX]{x1D6FB}_{\Vert }^{2}u=p,\end{eqnarray}$$

where $B=Et^{3}/[12(1-\unicode[STIX]{x1D708}^{2})]$ is the bending energy (flexural rigidity) of the plate, $E$ is the material’s Young’s modulus and $\unicode[STIX]{x1D708}$ is the Poisson ratio of the material; $\unicode[STIX]{x1D6FB}_{\Vert }^{2}$ is the Laplacian operator in the $(x,z)$ coordinates (tangent to the base flow). Equation (3.2) models only bending of the plate, with stretching being assumed negligible; if stretching is judged to be significant, then the Föppl–von Kármán equations can be employed (Timoshenko & Woinowsky-Krieger Reference Timoshenko and Woinowsky-Krieger1959; Landau & Lifshitz Reference Landau and Lifshitz1986).

Using the dimensionless variables from (2.3), equation (3.2) becomes

(3.3) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}^{4}\frac{\unicode[STIX]{x2202}^{4}U}{\unicode[STIX]{x2202}X^{4}}+2\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D6FF}^{2}\frac{\unicode[STIX]{x2202}^{4}U}{\unicode[STIX]{x2202}X^{2}\unicode[STIX]{x2202}Z^{2}}+\unicode[STIX]{x1D716}^{4}\frac{\unicode[STIX]{x2202}^{4}U}{\unicode[STIX]{x2202}Z^{4}}=\frac{h_{0}^{4}\unicode[STIX]{x0394}p}{Bu_{\mathit{c}}}P.\end{eqnarray}$$

To obtain a dominant balance in (3.3), we must take $h_{0}^{4}\unicode[STIX]{x0394}p/(Bu_{\mathit{c}})=\unicode[STIX]{x1D6FF}^{4}$ or, specifically, this balance sets the characteristic deformation scale to be

(3.4) $$\begin{eqnarray}u_{\mathit{c}}=\frac{h_{0}^{4}\unicode[STIX]{x0394}p}{B\unicode[STIX]{x1D6FF}^{4}}=\frac{w^{4}\unicode[STIX]{x0394}p}{B}.\end{eqnarray}$$

Hence, $\unicode[STIX]{x1D6FD}\equiv u_{\mathit{c}}/h_{0}=h_{0}^{3}\unicode[STIX]{x0394}p/(B\unicode[STIX]{x1D6FF}^{4})$ ; typical values of $\unicode[STIX]{x1D6FD}$ are estimated in table 1. With $u_{\mathit{c}}$ determined, equation (3.3) gives, to the leading order,

(3.5) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}^{4}U}{\unicode[STIX]{x2202}X^{4}}=P(Z).\end{eqnarray}$$

Physically, this asymptotic limit can be interpreted as follows: since $w\ll \ell$ , the variation of the height of the channel in the streamwise direction occurs over a much longer length scale than the variation of the height in any cross-section. In other words, the deflection of any infinitesimal slice (perpendicular to the flow-wise axial direction) of the top wall does not affect infinitesimal slices nearby. Thus, the Kirchhoff–Love plate theory reduces to the Euler–Bernoulli beam theory for each infinitesimal spanwise slice of the top wall. This limit should be contrasted to the case when $\unicode[STIX]{x1D6FF}\sim \unicode[STIX]{x1D716}$ (and the elastic top wall is not clamped along $X=\pm 1/2$ ) in which case bending (and, potentially, tension) in the flow-wise direction dominate (Gomez et al. Reference Gomez, Moulton and Vella2017).

Table 1.

Values of the physical parameters for the flow-rate-controlled PDMS microchannels experimental system S4 of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013), which most closely matches the assumptions of our theoretical developments.

Of course, at this order in the asymptotic analysis, the beam equation (3.5) cannot satisfy the boundary conditions along $Z=0$ or $Z=1$ (the channel’s entrance and exit planes). In experiments, the plate is also clamped along $Z=0$ and $Z=1$ . From (3.3), it is evident that a boundary-layer calculation can be done by the rescaling $Z\mapsto (\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6FF})Z$ , which keeps all terms on the left-hand side of (3.3). Thus, any corrections due to clamping at $Z=0$ and $Z=1$ are localized to layers of width $O(\unicode[STIX]{x1D716}/\unicode[STIX]{x1D6FF})$ , which are a small correction to the system’s response in the central part of the channel, given our assumed ordering of small parameters (2.6).

We take the plate to be clamped at $X=\pm 1/2$ , so (3.5) is subject to the boundary conditions

(3.6a,b ) $$\begin{eqnarray}U(\pm 1/2,Z)=0,\quad \left.\frac{\unicode[STIX]{x2202}U}{\unicode[STIX]{x2202}X}\right|_{X=\pm 1/2}=0.\end{eqnarray}$$

Integrating (3.5) four times with respect to $X$ and enforcing the boundary conditions from (3.6a,b ), we find that

(3.7) $$\begin{eqnarray}U(X,Z)={\textstyle \frac{1}{24}}P(Z)\left(X+{\textstyle \frac{1}{2}}\right)^{2}\left(X-{\textstyle \frac{1}{2}}\right)^{2}.\end{eqnarray}$$

Note that, because $u_{c}$ is set by bending via (3.4), we cannot take a limit as a parameter in (3.7) vanishes and recover the case of a rigid channel (i.e. $U=0$ ).

The top wall displacement given by (3.7) is a quartic polynomial of the cross-sectional coordinate $X$ , which corresponds to our assumption that the system is in a bending-dominated regime. Nevertheless, it should also be noted that a (different) quartic profile is consistent with recent experimental measurements of a similar model system (Ducloué et al. Reference Ducloué, Hazel, Thompson and Juel2017), in which bending and pre-tension can both be considered significant. In other literature (Gervais et al. Reference Gervais, El-Ali, Günther and Jensen2006; Cheung et al. Reference Cheung, Toda-Peters and Shen2012; Ozsun et al. Reference Ozsun, Yakhot and Ekinci2013; Raj & Sen Reference Raj and Sen2016), however, a parabolic (quadratic) profile has often been assumed, even without making a distinction between (or an evaluation of) bending-dominated versus stretching-dominated deformation.

It is often of interest to relate the cross-sectionally averaged displacement to the maximum displacement, since the latter is easier to measure experimentally. From (3.7) we find that

(3.8a,b ) $$\begin{eqnarray}\overline{U}(Z):=\int _{-1/2}^{+1/2}U(X,Z)\,\text{d}X=\frac{8}{15}U_{max}(Z),\quad U_{max}(Z):=\frac{1}{384}P(Z).\end{eqnarray}$$

In passing, we note that the prefactor in (3.8a ), when computed on the basis of the parabolic displacement approximation is $2/3$ rather than $8/15$ ; a 25 % difference. The predicted linear scaling of the maximum and cross-sectionally averaged displacements with the pressure is also consistent with estimates of the static deflection of plates by Timoshenko & Woinowsky-Krieger (Reference Timoshenko and Woinowsky-Krieger1959).

4 The flow rate–pressure drop relation

To find the flow rate through the deformed channel, we must evaluate the integral in (2.10) using the top wall shape stipulated in (3.1):

(4.1) $$\begin{eqnarray}Q=-\frac{1}{12}\frac{\text{d}P}{\text{d}Z}\int _{-1/2}^{+1/2}[1+\unicode[STIX]{x1D6FD}U(X,Z)]^{3}\,\text{d}X.\end{eqnarray}$$

This integral can be evaluated upon substituting the expression for $U$ from (3.7) to find

(4.2) $$\begin{eqnarray}Q=-\frac{1}{12}\frac{\text{d}P}{\text{d}Z}\left[1+\frac{\tilde{\unicode[STIX]{x1D6FD}}}{10}P(Z)+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{2}}{210}P(Z)^{2}+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{3}}{12\,012}P(Z)^{3}\right],\end{eqnarray}$$

where $\tilde{\unicode[STIX]{x1D6FD}}:=\unicode[STIX]{x1D6FD}/24=h_{0}^{3}\unicode[STIX]{x0394}p/(24B\unicode[STIX]{x1D6FF}^{4})$ . In passing, we note that (4.2) is a relationship of the form $Q=-\unicode[STIX]{x1D70E}(P)\,\text{d}P/\text{d}Z$ , for some function $\unicode[STIX]{x1D70E}$ that comes about from solving the elasticity problem, as discussed by Rubinow & Keller (Reference Rubinow and Keller1972).

For a constant imposed flow rate $Q$ , equation (4.2) is a separable first-order ordinary differential equation for the pressure distribution $P(Z)$ subject to $P(1)=0$ , whence

(4.3) $$\begin{eqnarray}Q=\frac{P(Z)}{12(1-Z)}\left[1+\frac{\tilde{\unicode[STIX]{x1D6FD}}}{20}P(Z)+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{2}}{630}P(Z)^{2}+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{3}}{48\,048}P(Z)^{3}\right].\end{eqnarray}$$

As $Q$ is constant, equation (4.3) is a polynomial, whose real positive root gives the dimensionless fluid pressure $P$ as a function of the streamwise coordinate $Z$ and the various dimensionless parameters. (Though it is possible to obtain the roots of the quartic polynomial in (4.3) by Ferrari’s procedure (see, e.g. Abramowitz & Stegun Reference Abramowitz and Stegun1972, § 3.8.3), the explicit formula for the sought positive real root is far too long to be presented here. It can, however, be found easily using Mathematica.) For $\tilde{\unicode[STIX]{x1D6FD}}=0$ (an ‘infinitely stiff’ plate), there is no deformation, and (4.3) reduces to the classical lubrication-theory result for a rectangular channel.

Since the last two terms in the bracket in (4.3) have fairly large denominators, they can be neglected for small as well as moderate values of $\tilde{\unicode[STIX]{x1D6FD}}P(Z)$ , the consequences of which are discussed in appendix A. Obviously, for larger values of $\tilde{\unicode[STIX]{x1D6FD}}P(Z)$ these higher-order terms become important.

4.1 Flow-rate-controlled versus pressure-drop-controlled problems

For a pressure-driven flow, $\unicode[STIX]{x0394}P=1$ under the non-dimensionalization $p=\unicode[STIX]{x0394}pP$ ; hence $P(0)=1$ . Then, the flow rate–pressure drop relation (4.3) can be evaluated at $Z=0$ to obtain

(4.4) $$\begin{eqnarray}12Q=1+\frac{\tilde{\unicode[STIX]{x1D6FD}}}{20}+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{2}}{630}+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{3}}{48\,048},\end{eqnarray}$$

which is a universal dimensionless relationship valid for any applied $\unicode[STIX]{x0394}p$ .

For a flow-rate-controlled situation, $Q=1$ under the non-dimensionalization $p=\unicode[STIX]{x1D707}\ell q/(h_{0}^{3}w)P$ ; hence $P(0)=\unicode[STIX]{x0394}P$ . Then, the flow rate–pressure drop relation (4.3) can be evaluated at $Z=0$ to obtain

(4.5) $$\begin{eqnarray}12=\unicode[STIX]{x0394}P\left[1+\frac{\tilde{\unicode[STIX]{x1D6FD}}}{20}\unicode[STIX]{x0394}P+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{2}}{630}(\unicode[STIX]{x0394}P)^{2}+\frac{\tilde{\unicode[STIX]{x1D6FD}}^{3}}{48\,048}(\unicode[STIX]{x0394}P)^{3}\right],\end{eqnarray}$$

which is a universal dimensionless relationship valid for any imposed $q$ .

Equations (4.4) and (4.5) can also be interpreted as consistency checks. For a given experiment (or numerical simulation) with an imposed (dimensionless) flow rate $Q$ or an imposed dimensionless pressure drop $\unicode[STIX]{x0394}P$ , the value of $\tilde{\unicode[STIX]{x1D6FD}}$ corresponding to the system must satisfy either the (4.4) or (4.5).

4.2 Comparison with previous models

Equation (4.3), which contains no fitting parameters, can be compared to equation (10) of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006), which incorporates a fitting parameter $\unicode[STIX]{x1D6FC}$ and has the dimensional form:

(4.6) $$\begin{eqnarray}q=\frac{h_{0}^{3}wp(z)}{12\unicode[STIX]{x1D707}(\ell -z)}\left[1+\frac{3}{2}\frac{\unicode[STIX]{x1D6FC}w}{Eh_{0}}p(z)+\left(\frac{\unicode[STIX]{x1D6FC}w}{Eh_{0}}\right)^{2}p(z)^{2}+\frac{1}{4}\left(\frac{\unicode[STIX]{x1D6FC}w}{Eh_{0}}\right)^{3}p(z)^{3}\right].\end{eqnarray}$$

Equation (4.6) was derived by Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006) under the assumption that the top wall is thick, and the strains within it, $u/w$ , are proportional to $p/E$ . The proportionality constant is $\unicode[STIX]{x1D6FC}$ , and it is unknown a priori. Although Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006) only discuss the case of a thick top wall that behaves like an elastic half-space, equation (4.6) has been used by others (see, e.g. Hardy et al. Reference Hardy, Uechi, Zhen and Kavehpour2009; Cheung et al. Reference Cheung, Toda-Peters and Shen2012; Raj et al. Reference Raj, DasGupta and Chakraborty2017) for microchannels with much thinner top walls that behave as plates. Consequently, when (4.6) is fit to flow pressure–drop data from microchannels with plate-like top walls, a dependence $\unicode[STIX]{x1D6FC}\propto (t/w)^{3}$ is observed (Raj et al. Reference Raj, DasGupta and Chakraborty2017).

For the bending-dominated analysis presented herein, we found (similarly to Gervais et al. Reference Gervais, El-Ali, Günther and Jensen2006) that $u/w\propto p/E$ , but on the basis of plate theory. Consequently, it is possible to compare our flow rate–pressure drop relation (4.3) to a generic relation such as (4.6) since they are both based on the top wall displacement scaling with the fluid pressure, i.e. $u/w\propto p/E$ . In order to discuss the similarities and differences between these two relations, we note that the dimensional form of (4.3) is

(4.7) $$\begin{eqnarray}\displaystyle q & = & \displaystyle \frac{h_{0}^{3}wp(z)}{12\unicode[STIX]{x1D707}(\ell -z)}\left[1+\frac{1}{480}\frac{w^{4}}{Bh_{0}}p(z)+\frac{1}{362\,880}\left(\frac{w^{4}}{Bh_{0}}\right)^{2}p(z)^{2}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\frac{1}{664\,215\,552}\left(\frac{w^{4}}{Bh_{0}}\right)^{3}p(z)^{3}\right].\end{eqnarray}$$

Clearly, the quantity in the brackets on the right-hand side of (4.6) can be compared with the quantity in the brackets on the right-hand side of (4.7). The caveat here is that we assume that (4.6) can be applied to the bending-dominated top wall deformation regime. In this comparison, then, the role of the parameter $\unicode[STIX]{x1D6FC}$ is to absorb the necessary details of the displacement profile and its scaling with various material and geometric constants. For example, the second terms in the brackets of (4.6) and (4.7) can be made to agree if we choose $\unicode[STIX]{x1D6FC}=(1/60)(w/t)^{3}(1-\unicode[STIX]{x1D708}^{2})$ . We cannot make all terms agree because (4.7) is based on the cross-sectional displacement profile $u(x,z)$ derived in § 3, while Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006) postulated that $u/w=\unicode[STIX]{x1D6FC}p/E$ .

Raj & Sen (Reference Raj and Sen2016) and Raj et al. (Reference Raj, DasGupta and Chakraborty2017) also observed that $\unicode[STIX]{x1D6FC}$ need not be a fitting parameter but, instead, must be related to the various geometric and elasticity constants involved in the problem. A difference between the approaches of Raj & Sen (Reference Raj and Sen2016) and Raj et al. (Reference Raj, DasGupta and Chakraborty2017) and ours is that their derivations are based on correlating the maximum of the cross-sectionally averaged displacement to the pressure via various approximations/assumptions. In particular, Raj & Sen (Reference Raj and Sen2016) make the assumption of a parabolic top wall deflection, which does not appear to be well substantiated by their measurements. Consequently, the (dimensional) ‘compliance parameters’, $f_{p}$ introduced by Raj & Sen (Reference Raj and Sen2016) and $\unicode[STIX]{x1D6FA}$ introduced by Raj et al. (Reference Raj, DasGupta and Chakraborty2017), are derived differently from (the dimensionless) $\tilde{\unicode[STIX]{x1D6FD}}$ introduced above on the basis of the perturbation expansion of the coupled fluid–structure interaction problem.

More importantly, however, Raj & Sen (Reference Raj and Sen2016) based their analysis on the shell-stretching correlation $u/w\propto (p/E)^{1/3}$ , which yields a model that cannot be compared to those based on linearly proportional displacement and pressure (such as those discussed above). On the other hand, the model flow rate–pressure drop relation of Raj et al. (Reference Raj, DasGupta and Chakraborty2017) was based on a plate-bending correlation, so it can be directly compared to (4.7). Unfortunately, there are misprints in the derivation of equation (12) of Raj et al. (Reference Raj, DasGupta and Chakraborty2017), i.e. the model flow rate–pressure drop relation therein. Specifically, the compliance parameter $\unicode[STIX]{x1D6FA}$ as given in equation (10) of Raj et al. (Reference Raj, DasGupta and Chakraborty2017) should read $\unicode[STIX]{x1D6FA}=(1/99)E^{-1}(w/h_{0})(w/t)^{3}(1-\unicode[STIX]{x1D708}^{2})$ and the flow rate–pressure drop relation from equation (12) of Raj et al. (Reference Raj, DasGupta and Chakraborty2017) should read

(4.8) $$\begin{eqnarray}q=\frac{h_{0}^{3}wp(z)}{12\unicode[STIX]{x1D707}(\ell -z)}\left[1+\frac{3}{2}\unicode[STIX]{x1D6FA}p(z)+\unicode[STIX]{x1D6FA}^{2}p(z)^{2}+\frac{1}{4}\unicode[STIX]{x1D6FA}^{3}p(z)^{3}\right].\end{eqnarray}$$

On comparing (4.8) to (4.6), it follows that $\unicode[STIX]{x1D6FC}=(1/99)(w/t)^{3}(1-\unicode[STIX]{x1D708}^{2})$ . Unsurprisingly, this comparison yields a constant $\unicode[STIX]{x1D6FC}$ , which is no longer a fitting parameter. The latter prediction differs from ours above in the prefactor $(1/99)\approx 0.011$ versus $(1/60)\approx 0.017$ . Furthermore, as was the case with (4.6), equation (4.8) cannot be made to agree with (4.7) because (4.8) was derived following the approach of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006), while (4.7) was derived using the actual cross-sectional displacement profile $u(x,z)$ .

5 Illustrated examples and comparison with experiments

First, consider the leading-order axial velocity component (2.8). Clearly, because $P=P(Z)$ and $H=H(X,Z)$ , the deformation of the channel’s top wall introduces dependence upon both of the coordinates (i.e. $X$ and $Z$ ) that are not present in the leading-order velocity in a rigid channel, namely $V_{Z}(X,Y,Z)=6Q(1-Y)Y$ . The effect of the top wall’s elasticity on the axial velocity is illustrated in figure 2. Specifically, it is evident that the cross-sectional shape’s deformation introduces variability in the $X$ -direction, which leads to faster flow near the channel’s centreline (compared to the rest of the cross-section), consistent with the numerical simulations shown in figure 3 of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006).

Figure 2.

Contours of the axial velocity, $V_{Z}$ , as a function of the $(X,Y)$ coordinates in the $Z=0.25$ cross-section of the channel with $Q=1$ and for different $\tilde{\unicode[STIX]{x1D6FD}}$ values, computed from (2.8) with the pressure given by the appropriate real root of (4.3) and the channel height given by (3.1) and (3.7). (a) Solid top wall ( $\tilde{\unicode[STIX]{x1D6FD}}=0$ ). (b) Elastic top wall with $\tilde{\unicode[STIX]{x1D6FD}}=0.5$ . (c) Elastic top wall with $\tilde{\unicode[STIX]{x1D6FD}}=1$ .

Second,the channel’s top wall deformation is illustrated in figure 3. Indimensionless variables, figure 3(a) shows a surface plot of the top wall displacement as a function of the spatial coordinates, while figure 3(b) shows the maximum displacement of the top wall as a function of the flow rate $Q$ . The shapes of the latter curves agree qualitatively with those in figures 5 and 7 of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006), respectively.

Figure 3.

(a) $H(X,Z)$ computed from (3.7) with $P(Z)$ computed from (4.3), for $\tilde{\unicode[STIX]{x1D6FD}}=2$ and $Q=1$ ; thick dashed line represents $1+\unicode[STIX]{x1D6FD}U_{max}(Z)$ , thin dotted lines are a guide to the eye. (b) Maximum displacement, i.e. $H(0,0)-1=\unicode[STIX]{x1D6FD}U(0,0)$ calculated from (3.7) with $P(Z)$ computed from (4.3), as function of the flow rate $Q$ and different values of $\tilde{\unicode[STIX]{x1D6FD}}$ .

Third, the nonlinear flow rate–pressure drop curves for a deformable channel are illustrated in figure 4(a) for different values of $\tilde{\unicode[STIX]{x1D6FD}}$ . Their shapes agree qualitatively with Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006) (figure 9 therein). From figure 4(a), it is clear that, for $\tilde{\unicode[STIX]{x1D6FD}}\neq 0$ (a deformable channel), the value of $Q$ for a given $\unicode[STIX]{x0394}P$ can be significantly larger as $\unicode[STIX]{x1D6FD}$ is increased, i.e. for softer channels. Meanwhile, figure 4(b) shows the pressure in the microchannel as a function of the axial coordinate using the implicit expression (4.3). Clearly, a contribution of flow–elasticity coupling is to make the pressure distribution in the channel nonlinear. Furthermore, since the channel cross-sectional area increases due to the top wall’s deformation, a smaller pressure gradient is needed to achieve the same volumetric flow rate.

Figure 4.

(a) Pressure drop across the channel, $\unicode[STIX]{x0394}P=P(0)$ , computed via (4.3), as a function of the flow rate $Q$ for different values of $\unicode[STIX]{x1D6FD}$ . (b) Pressure $P$ as a function of the axial coordinate $Z$ computed by inverting (4.3);  $Q=1$ .

Next, we present comparisons between the predictions of our theoretical analysis and experimental results previously published in the literature. To this end, table 1 summarizes the values of the physical parameters in the flow-rate-controlled PDMS microchannels experiment, which best matches the assumptions of our theoretical developments, labelled ‘S4’ in Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013). The ranges of $\unicode[STIX]{x1D6FD}$ and $\tilde{\unicode[STIX]{x1D6FD}}$ , which were computed using (4.4), correspond to the range of flow rates considered by Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013), namely $q=10{-}50~\text{ml}~\text{min}^{-1}$ , or $Q=0.1{-}0.15$ . From table 1, which also shows the values of the dimensionless quantities that we have introduced, we see that $\unicode[STIX]{x1D6FF}\approx 0.144$ while $\unicode[STIX]{x1D716}\approx 0.02$ , which satisfies our assumed asymptotic ordering: $0<\unicode[STIX]{x1D716}\ll \unicode[STIX]{x1D6FF}\ll 1$ . Regarding our assumption that plate-bending theory captures the mechanics of the top wall deformation, we note, from table 1, that $t/w=0.11\ll 1$ , while $\max _{x,z}|u(x,z)|/t\leqslant 0.43<1$ based on the results of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013) (with $0.43$ being the value at the largest flow rate considered therein).

To make a quantitative comparison between theory and experiment, figure 5 shows a plot of the dimensional flow rate–pressure drop relation (4.7) based on our asymptotic theory, the corresponding experimental measurements from the S4 system of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013) and a reference line corresponding to the lubrication-theory result for a solid rectangular channel. Clearly, there is good agreement between theory and experiment, especially with respect to the shape of the flow rate–pressure drop curve. In particular, the cross-over from a linear (at low $q$ ) to a nonlinear (at high $q$ ) regime can be observed in the plot.

Figure 5.

Dimensional flow rate–pressure drop relationship: comparison between our asymptotic theory (solid curve) (i.e. (4.7) evaluated at $z=0$ ), including a shaded trust region based on taking $E=1.6\times 10^{6}\pm 25\,\%$ MPa, and the experimental data from the S4 system of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013) (symbols) with error bars smaller than the symbols according to Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013). The dotted line represents the standard lubrication-theory linear flow rate–pressure drop relation for a rigid rectangular channel, i.e. $\unicode[STIX]{x0394}p=12\unicode[STIX]{x1D707}\ell q/(h_{0}^{3}w)$ .

The theoretical curve shown in figure 5 systematically underpredicts the pressure drop by approximately 5 %. We attribute this small discrepancy to two factors. First, the deformation of the S4 microchannel is not insignificant: as reported in table 1 of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013), $\max _{x,z}|u(x,z)|/t\approx 0.43$ . Such a maximum displacement, given the top wall’s thickness, might be considered at the edge of applicability of the Kirchhoff–Love (linearly elastic) plate theory. Hence, the agreement might be improved if a different top wall elasticity model, which accounts for the non-trivial wall thickness, were used. Second, the elasticity parameters of the PDMS microchannels are not well characterized by Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013). We have estimated a value of $E\approx 1.6~\text{MPa}$ using their description of the experimental procedure, cross-referencing against the data on PDMS material properties tabulated by Johnston et al. (Reference Johnston, McCluskey, Tan and Tracey2014) and treating the hydrostatic loading data (Ozsun et al. Reference Ozsun, Yakhot and Ekinci2013, figure 1(e,f)) as a ‘bulge test’ (Small & Nix Reference Small and Nix1992). However, there is inherent uncertainty in such an estimate, which is why we have added a confidence region (corresponding to a $\pm 25\,\%$ variation in the Young’s modulus $E$ ) around the theoretical curve in figure 5. This uncertainty appears to capture most of the (small) discrepancy between theory and experiment. Therefore, it is clear that figure 5 shows that a quantitative prediction of the flow rate–pressure drop curve for a long shallow microchannel is possible without fitting parameters.

Finally, based on the discussion § 4.2, it might be appropriate to think of the fitting parameter $\unicode[STIX]{x1D6FC}$ of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006) as constant given by $\unicode[STIX]{x1D6FC}=(1/60)(w/t)^{3}(1-\unicode[STIX]{x1D708}^{2})$ . Then, for the S4 experiments of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013), we find that fitting the data to (4.6) evaluated at $z=0$ yields $\unicode[STIX]{x1D6FC}\approx 5.91$ . (Here, we have used Mathematica’s built-in nonlinear least-squares subroutine FindFit. Note, however, that the best-fit value is sensitive to small errors introduced in digitizing the data from Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013).) Meanwhile, $\unicode[STIX]{x1D6FC}=(1/60)(w/t)^{3}(1-\unicode[STIX]{x1D708}^{2})\approx 7.68$ using the values given in table 1.

6 Discussion and conclusions

Using asymptotic techniques, we have derived a flow rate–pressure drop relation for low Reynolds number flow in a shallow deformable channel. Our result, which contains no fitting parameters, can be applied to the design and construction of soft microfluidic devices, where the flow and the elastic deformation of the domain are coupled. Specifically, we treated the case of bending-dominated deformation, which is governed by the plate equation (3.2). This restricts our analysis to microchannels for which the maximum deformation $u_{max}:=\max _{x,z}|u(x,z)|$ is much smaller than the thickness of the top wall, which is, in turn, much smaller than the channel’s width: $u_{max}\ll t\ll w$ .

To make clear the limitations of our model, consider first the case when $t\gg w$ , i.e. the top wall is much thicker than the width of the channel, as is the case in the experiments of Gervais et al. (Reference Gervais, El-Ali, Günther and Jensen2006). Then, the top wall behaves like an elastic half-space rather than a plate. An extension of the present work would require us to show a similar decoupling between the flow-wise and transverse deformation, as was done in § 3, perhaps adapting the solution for a uniform load distributed over a finite width of an elastic half-space (Johnson Reference Johnson1985, § 2.4). Now, consider the case when $u_{max}\gg t$ , i.e. the top wall’s deformation is much larger than its thickness, as is the case in some of the experiments of Ozsun et al. (Reference Ozsun, Yakhot and Ekinci2013). Now, the top wall behaves like a membrane rather than a plate. Once again, an extension of the present work would be to revisit § 3 starting from a model of stretching-dominated deformation, perhaps adapting the approximate solution from Timoshenko & Woinowsky-Krieger (Reference Timoshenko and Woinowsky-Krieger1959, Art. 101). Each case would require a careful re-calculation of the flow rate–pressure drop relation and appropriate comparisons against applicable experiments. Unfortunately, at this time, few reliable measurements of the top wall displacement of microchannels under hydrodynamic conditions have been published.

Returning to our results for the bending-dominated case, in dimensional variables and keeping only the leading-order elasticity contribution in (4.7), we can summarize the key result:

(6.1) $$\begin{eqnarray}q\approx \frac{h_{0}^{3}w\unicode[STIX]{x0394}p}{12\unicode[STIX]{x1D707}\ell }\left[1+\frac{3}{160}\left(\frac{w}{t}\right)^{3}\left(\frac{w}{h_{0}}\right)\left(\frac{\unicode[STIX]{x0394}p}{E}\right)\right],\end{eqnarray}$$

where $q$ is the flow rate, $\unicode[STIX]{x0394}p$ is the pressure drop, $w$ is the channel width, $h_{0}$ is the undeformed channel height, $\ell$ is the channel length, $t$ is the top wall’s thickness and $E$ is the Young’s modulus of the material (assuming an incompressible material, i.e. a Poisson ratio $\unicode[STIX]{x1D708}=1/2$ ). In particular, equation (6.1) highlights the strong dependence on the plate geometry through the factor of $(w/t)^{3}$ , which is assumed $\gg 1$ in our plate-bending model. Note that the ratio outside the square brackets in (6.1) represents the leading-order term (in $h_{0}/w\equiv \unicode[STIX]{x1D6FF}$ ) in the flow rate–pressure drop relation for a shallow rectangular channel.

However, the boundary conditions in the transverse $x$ -direction are not enforced in our lubrication-theory-based approach. This defect comes from the shallowness assumption that allows us to neglect $O(\unicode[STIX]{x1D6FF}^{2})$ terms in the asymptotic expansion. A potential remedy is to relax the assumption that $\unicode[STIX]{x1D6FF}\ll 1$ and solve the full two-dimensional leading-order flow problem using the domain perturbation technique (Van Dyke Reference Van Dyke1975; Lebovitz Reference Lebovitz1982). An alternative approach, consistent with lubrication theory, is presented in appendix B. Finally, future work might include computing the remaining two components of the leading-order velocity field. This three-dimensional case requires further consideration because, as shown by Lauga, Stroock & Stone (Reference Lauga, Stroock and Stone2004), if the cross-section is varying in the flow-wise direction, then the flow cannot be planar, so the continuity equation alone does not yield the remaining velocity component as in standard lubrication theory. It would also be of interest to compute higher-order perturbative corrections to the lubrication theory velocity profile given in § 2 following, e.g. Tavakol et al. (Reference Tavakol, Froehlicher, Holmes and Stone2017).