Soft hydraulics, which addresses the interaction between an internal flow and a compliant conduit, is a central problem in microfluidics. We analyze Newtonian fluid flow in a rectangular duct with a soft top wall at steady state. The resulting fluid–structure interaction is formulated for both vanishing and finite flow inertia. At the leading-order in the small aspect ratio, the lubrication approximation implies that the pressure only varies in the streamwise direction. Meanwhile, the compliant wall's slenderness makes the fluid–solid interface behave like a Winkler foundation, with the displacement fully determined by the local pressure. Coupling flow and deformation and averaging across the cross section leads to a one-dimensional reduced model. In the case of vanishing flow inertia, an effective deformed channel height is defined rigorously to eliminate the spanwise dependence of the deformation. It is shown that a previously used averaged height concept is an acceptable approximation. From the one-dimensional model, a friction factor and the corresponding Poiseuille number are derived. Unlike the rigid duct case, the Poiseuille number for a compliant duct is not constant but varies in the streamwise direction. Compliance can increase the Poiseuille number by a factor of up to four. The model for finite flow inertia is obtained by assuming a parabolic vertical variation of the streamwise velocity. To satisfy the displacement constraints along the edges of the channel, weak tension is introduced in the streamwise direction to regularize the Winkler-foundation-like model. Matched asymptotic solutions of the regularized model are derived.

## I. INTRODUCTION

As Frank M. White notes in Sec. 3–3.3 of his iconic *Viscous Fluid Flow* textbook,^{2} “fully developed duct flow is equivalent to a classic Dirichlet problem, [thus] it is not surprising that an enormous number of exact solutions are known.” He (and, also, Bruus^{3}) summarized an elegant set of solutions for such unidirectional flows, which are *exact* solutions of the incompressible Navier–Stokes equations at *any* Reynolds number, including “limaçon-shaped ducts, for example, [which] are not commercially available at present,”^{2} adding a bit of humor to this topic. The general result is that, for all such duct flows, the volumetric flow rate *q* is related to the axial pressure gradient, $\u2212dp/dz$, via some (likely complicated) function of the cross-sectional geometry. This result is the cornerstone of *hydraulics*, i.e., “the conveyance of liquids through pipes and channels” (per Oxford Languages—the provider of Google's English dictionary), a topic that is taught to undergraduate students.^{4}

Now, however, what if the duct were manufactured from a soft material so that the local hydrodynamic pressure changes the cross-sectional area? Such problems have a time-honored history in biomechanics^{5–7} but not so much in hydraulics. Nevertheless, with the emergence of *microfluidics,*^{3,8–11} the hydraulics of compliant ducts manufactured from polymeric materials has become a central problem at the intersection of fluid mechanics and soft matter physics. To develop a theory of *soft hydraulics*, we must understand steady fluid–structure interactions (FSIs). FSIs between external or internal flows (either viscous or inviscid) and elastic structures, as well as the linear stability of such coupled mechanics problems, is also a well-developed research subject,^{12} including fast progress in the last decade.^{13} While FSI topics such as aeroelasticity^{14} and moderate-Reynolds-number blood flow in large arteries^{6} are now quite classical, the mechanical interaction between *slow* viscous flows and *compliant* conduits^{15} has opened new avenues of FSI research,^{13,16} both at the microscale, e.g., lab-on-a-chip applications,^{15,17} and at the macroscale, e.g., soft robotics applications.^{18,19}

To this end, in this paper, soft hydraulics and its mathematical formulation are first introduced in Sec. II. Then, the discussion bifurcates into the case of a vanishing Reynolds number (Sec. III) and the case of a finite Reynolds number (Sec. IV). We review the key recent results regarding flows in compliant ducts of initially rectangular cross section. Then, within each of Secs. III and IV, we show how to consistently reduce these inherently three-dimensional (3D) problems to two-dimensional (2D) problems^{20} and, eventually, to one-dimensional (1D) models (that only involve axial, or streamwise, variations).^{21} Specifically, in Sec. III, we compare our consistent formulation with previous spanwise-averaged (i.e., over *x*, see Fig. 1) models and ascertain the accuracy of the previous approach to the hydraulic predictions. Toward this end, in Sec. III D, we introduce a generalization of the laminar flow friction factor suitable for quantifying the effect of compliance in soft hydraulic systems. In Sec. IV, we address the issue of Reynolds number dependence (flow acceleration), which is a novel contribution of our work to the field of soft hydraulics. However, the model breaks down beyond a certain Reynolds number, requiring a regularization (Sec. IV B), which leads to an interesting singular perturbation problem (solved in Appendix A). Finally, conclusions and avenues for future work are discussed in Sec. V.

## II. PRELIMINARIES, NOTATION, AND PROBLEM STATEMENT

Consider a soft-walled microchannel (initially a rectangular duct), which exhibits flow-induced deformation due to Newtonian fluid flow through it.^{22} Denote the channel's undeformed height, width, length, and top wall thickness by *h*_{0}, *w*, $\u2113$, and *t*, respectively, as in Fig. 1. Further, introducing the aspect ratios, $\delta =h0/w$ and $\u03f5=h0/\u2113$, we say that the microchannel is long and shallow^{23} if $\u03f5\u226a\delta \u226a1$. This kind of compliant duct is a common outcome of rapid microfluidic device fabrication via soft lithography.^{24,25} In the following analysis, the top wall's deformation is dominant, and thus it is the only deformation of interest.^{23,26} Denote the deformed cross-sectional height by *h*(*x*, *z*) and assume the smallness of the aspect ratios still holds in the deformed microchannel, i.e., the deformed channel height is such that $h(x,z)\u226aw\u226a\u2113$. Note that $h(x,z)=h0+uy(x,z)$ is the deformed channel height, where $uy(x,z)$ is the displacement of the fluid–solid interface.

The incompressible Navier–Stokes (iNS) equations at steady state govern the flow within the duct. The velocity field is denoted $v=(vx,vy,vz)$ in Cartesian coordinates. To make iNS dimensionless, let us introduce the following dimensionless variables^{23,26} (denoted by capital letters),

The characteristic velocity and pressure scales $Vc$ and $Pc$, respectively, are discussed below. Under this nondimensionalization, the leading-order terms (in *ϵ*) left in iNS are^{1,26}

The fluid domain is defined as the deformed conduit: ${(X,Y,Z)|\u22121/2<X<+1/2,0<Y<H(X,Z),0<Z<1}$, in terms of the dimensionless variables. Here, $Re\u0302$ is the modified Reynolds number, defined as $Re\u0302=\u03f5Re=\u03f5\rho Vch0/\mu $, where *ρ* and *μ* are the fluid's density and dynamic viscosity, respectively. Equation (5) relates the characteristic pressure and velocity scales as $Pc=\mu \u2113Vc/h02$. Equations (3) and (4) indicate that, at the leading order in *ϵ* and *δ*, the hydrodynamic pressure *P* is only a function of the streamwise location *Z*, as in classical hydraulics problems.^{2} Importantly, however, in this soft hydraulics problem, the hydraulic resistance [set by the cross-sectional shape and area^{2,3} via Eq. (5)] is *not* constant and also varies with *Z*.

The following discussion begins with the case of $Re\u0302\u21920$ (in Sec. III), i.e., flow with negligible inertia. In this case, we consider two different mechanical responses of the compliant microchannel's wall, for which analytical solutions, based on the notion of a slowly varying^{27} unidirectional flow solution,^{28} are available in the literature.^{23,26,29–32} For both types of mechanical response, the previous solutions yield a 3D model, in which the axial flow profile $VZ=VZ(X,Y,Z)$ and the top wall shape $H=H(X,Z)$ are coupled via the hydrodynamic pressure *P*(*Z*). Our goal here is to first construct and validate reduced 2D models by “removing” the *X* dependence in a suitably rigorous way, so that $H=H(Z)$ only. Upon accomplishing this reduction, averaging the 2D model over *Y* yields a 1D model in which $H=H(Z)$ and $P=P(Z)$ are the remaining dependent variables. Therefore, when we extend the model to account for $Re\u0302=O(1)$ (in Sec. IV), i.e., to flow with moderate inertia, it suffices to consider just one reduced model (instead of each mechanical response individually).

## III. NEGLIGIBLE FLOW INERTIA: $Re\u0302\u21920$

### A. Effective deformed channel height

Neglecting the inertia of the flow by taking $Re\u0302\u21920$ in Eq. (5), we find that the axial velocity *V _{Z}*, subject to the no-slip boundary condition (BC) at the walls, has a parabolic variation along the height of the duct (

*Y*-direction),

At steady state, the flow rate is

Equation (8) can satisfy either one or two pressure boundary conditions (BCs). On the one hand, if the flow rate is controlled, then we can enforce $Q=q/q=1$ [i.e., take $Pc=\mu q\u2113/(wh03)$ in the nondimensionalization] and set the outlet pressure to gauge, i.e., $P(Z=1)=0$. On the other hand, if the pressure drop $\Delta P=P(Z=0)\u2212P(Z=1)$ is controlled, then enforcing $P(0)=p(0)/\Delta p=1$ (i.e., taking $Pc=\Delta p$ in the nondimensionalization) is now also a BC, in addition to $P(Z=1)=0$, from which *Q* is determined like an eigenvalue. Thus, in principle, the dimensionless flow rate *Q* and the dimensionless pressure drop $\Delta P$ are not independent,^{23} and we do not specify the flow regime *a priori* to make our results general. Either way, the pressure distribution in the duct can be determined by integrating Eq. (8) in *Z*, as long as the shape of fluid–solid interface, i.e., *H*(*X*, *Z*) is known.

Before we introduce expressions for *H*(*X*, *Z*), recall that, in a wide rigid rectangular duct, the relation between the pressure gradient and the flow rate is set by a Poiseuille-like law,^{2}

Thus, for a clearer comparison, it is helpful to transform Eq. (8) back into the dimensional form as

In order to consistently rewrite Eq. (10) in the form of a Poiseuille-like law (9), we define the effective channel height as

Then, Eq. (10) can be rewritten as

Note that the corresponding dimensionless effective channel height is

Equation (12) can be viewed as a generalization of the Poiseuille-like law (for a wide rigid rectangular duct) to a variable-height microchannel. From another perspective, using the axially varying height $he(z)$ in Eq. (12) eliminates the spanwise *x*-dependence of *h*(*x*, *z*). Then, since the velocity was already averaged across the cross section (to introduce *q*), the original 3D model has been reduced to an effective 1D model. Note that *h _{e}* is meaningful only when speaking of the relation between

*q*and $dp/dz$, both of which only vary with

*z*. This fact does

*not*mean that the velocity field is also 1D (it still depends on

*both y*and

*z*, thus remaining 2D). The effective height concept will be used to evaluate the accuracy of previous empirically motivated reduced-order models.

In particular, in the original studies using 1D models, such as those proposed by Gervais *et al.*^{22} and Hardy *et al.*,^{33} the *average* deformed channel height,

is used in Eq. (12) instead of $he(z)$. The corresponding dimensionless averaged channel height is

It should be clear, however, that $h\xaf$ (or $H\xaf$) from Eq. (14) [or Eq. (15)] is *not* equal to *h _{e}* (or

*H*) from Eq. (11) [or Eq. (13)]. Importantly, the averaging approach (introducing $h\xaf$ instead of

_{e}*h*) leads to an inconsistency in the reduced model because if we replace $he3(z)$ with $h\xaf3(z)$ in Eq. (12), then it is no longer equivalent to Eq. (10), which was rigorously derived by integrating the leading-order iNS (2)–(5). In the present work, our goal is to determine how this inconsistency affects the hydraulic predictions.

_{e}To finish the derivation, we must specify *H*(*X*, *Z*). In the present context of soft hydraulics, *H*(*X*, *Z*) is determined by solving an appropriate solid mechanics (elasticity) problem. Our assumption is that the deformed microchannel remains long and shallow [$h(x,z)\u226aw\u226a\u2113$] so that $uy\u226aw$ [recall that $h0\u226aw$ and $h(x,z)=h0+uy(x,z)$]. If the top wall is thick enough, with $w\u2272t\u226a\u2113$, then the shallowness and slenderness of the deformed channel enforces small-strain deformation and allows the use of linear elasticity. If the top wall is thin with $t\u2272w\u226a\u2113$, we require that $maxx,zuy\u226at$ to make the linear elastic theory applicable.^{23,29,30} However, regardless of the wall thickness, as long as $t\u226a\u2113$, the original 3D elasticity problem can be reduced to a 2D one. Here, we only briefly outline the reasons for the statement, and the reader is directed to Ref. 26 for the detailed analysis. First, the lubrication approximation implies^{2} that the shear stress $\tau yz\u226ap$. Since the tractions are continuous across the fluid–solid interface, we can thus infer that $\sigma yz\u226a\sigma yy$ where *σ _{yz}* and

*σ*are components of the Cauchy stress. Then, by examining the momentum balance in the solid, it is concluded that the dominant components of stress are in the cross-sectional (

_{yy}*x*,

*y*) plane. As a consequence, the deformation profiles at different streamwise (

*z*-locations) decouple from each other, leading to a local deformation–pressure relation.

Now, from Eqs. (3) and (4), $P=P(Z)$ only; thus, *P* acts uniformly at each axial (*X*, *Y*) cross section to deform the top wall. Therefore, the spanwise deformation is determined by the local pressure *P*(*Z*), and we can express the deformed duct shape as

where $\lambda :=Uc/h0$, with $Uc$ being the characteristic displacement of the top wall, is a dimensionless group that captures the compliance of the top wall. Restating the above-mentioned slenderness assumptions, we must require that $\lambda \u226a1/\delta $ for lubrication theory and linear elasticity to be applicable. The spanwise profile *F*(*X*) is obtained by solving the corresponding elasticity problem in the (*X*, *Y*) cross section of the duct.^{23,26,29,30} Also, note that Eq. (16) is *not* an assumption but a consequence of the asymptotic reduction of the elasticity problem for a long and slender microchannel. Since the analysis (summarized above) only involves balancing the momentum equation in the solid, it holds for any boundary conditions. However, the boundary conditions play a role in determining the actual deformation field, leading to different expressions for *F*(*X*).

Note that Eq. (16) takes the form of the deformation of a soft interface on a Winkler foundation,^{34,35} but now the foundation's (dimensionless) “spring stiffness” is given by $\lambda F(X)$. Winkler-foundation-like relations between pressure and deformation arise in a number of soft lubrication problems,^{36,37} including particles near elastic substrates,^{38,39} slider bearings,^{40} and rollers.^{41} The analogy becomes even stronger upon introducing the concept of averaged deformed channel height. Specifically, substituting Eq. (16) into Eqs. (13) and (15), respectively, we obtain explicit expressions for $He(Z)$ and $H\xaf(Z)$ as

and

Then, from Eq. (18), the now-constant (dimensionless) spring stiffness in the analogy to a Winkler foundation is $\xi =I1\lambda $. Here, the coefficients $Ii$ are defined as

Interestingly, observe that $H\xaf$ in Eq. (18) is simply the one-term Taylor-series approximation of *H _{e}* from Eq. (17) in terms of $\lambda \u226a1$. However, our analysis does not require $\lambda \u226a1$; in fact, $\lambda =O(1)$ is possible. Linear elasticity only requires that $\lambda \u226a1/\delta $ (as discussed by Wang and Christov

^{26}and Shidhore and Christov

^{29}) Thus, we would like to determine if the approximation in going from Eq. (17) to Eq. (18) is a valid one.

### B. Flow rate–pressure drop relation

To obtain the general form of the flow rate–pressure drop relation in a soft hydraulic conduit, we return to the dimensionless form of Eq. (8), namely,

Since $Q=const.$ in steady flow, upon substituting Eq. (17) into Eq. (20), we obtain a separable first-order ordinary differential equation (ODE) for *P*(*Z*). The solution, subject to $P(1)=0$, is

As discussed in Sec. III A, previous empirical studies used $H\xaf$ in place of *H _{e}*. In this case, substituting Eq. (18) into Eq. (20), and solving the corresponding ODE, yields an explicit expression for the pressure distribution,

As mentioned in Sec. III A, we may either consider a flow-controlled situation, in which *Q *=* *1 and $\Delta P=P(0)$ is found implicitly from Eq. (21) or explicitly from Eq. (22). Meanwhile, in the pressure-controlled regime, we enforce $P(0)=1$ and compute *Q* directly,

Equation (22) is essentially the same model derived by Gervais *et al.*^{22} However, in said work, $I1\lambda $ was taken as an unknown parameter, denoted as *α*, which was calibrated against experiments. However, our Eq. (22) is parameter-free because both *λ* and $I1$ are known from solving a suitable elasticity problem. Therefore, our approach eliminates the ambiguity, pointed out by Hardy *et al.*,^{33} of what unknown dependencies “hide” in *α*.

Note, however, that even if Eq. (18) is the one-term Taylor-series approximation to Eq. (17), this is not true for the flow rate–pressure drop relations in Eqs. (22) and (21), respectively. Therefore, we must determine how well *P*(*Z*) based on the averaged channel height approximates *P*(*Z*) based on the effective channel height. It is reasonable to conjecture that, due to the restriction to small strains required by linear elasticity, the two expressions should be in close agreement. To substantiate this conjecture, we proceed to quantify the difference between Eqs. (21) and (22) to obtain insights into the error committed in the formulation based on the averaged channel height. To this end, we apply the methodology established in this subsection to two types of common microchannel wall deformations considered in the literature: a microchannel with a thick top wall (Sec. III C 1) and a microchannel with a thinner, plate-like top wall (Sec. III C 2).

### C. Illustrated examples

#### 1. Duct with thick compliant top wall

First, we analyze the case of an initially rectangular duct with three compliant walls embedded in a thick soft structure. The channel's shallowness makes the deformation of the sidewall negligible compared with that of the top wall. Thus, the schematic diagram in Fig. 1 still applies. The corresponding steady 3D FSI problem was solved by Wang and Christov.^{26} To summarize their key conclusions: although a solution was obtained for any *t*/*w*, it was shown that, for $t/w\u22731.5$, the “thick” limit ($t2/w2\u226b1$) is achieved, and a simple analytical Fourier series solution can be written down for the deformed channel's top wall,

where we have defined $Am:=2m\pi [1\u2212(\u22121)m]$ and $E\xaf:=E/(1\u2212\nu 2)$, with *E* being Young's modulus and *ν* the Poisson's ratio.

From Eq. (25), we can determine the function $F(X)\u2261F(x/w)=f(x)$ introduced in Eq. (16). The corresponding values of $Ii$, defined in Eq. (19), are computed and summarized in Table I.

The compliance parameter *λ* emerges naturally from the nondimensionalization of Eq. (24),

Substituting *λ* and $Ii$ into Eqs. (21) and (22), respectively, we are ready to make a comparison between the two formulations. We observe that the pressure distribution depends nonlinearly upon *λ*, as illustrated in Fig. 2. The total pressure drop $\Delta P=P(0)$ decreases with *λ*, and a strong pressure gradient develops near the outlet. Notably, even with *λ* varying by three orders, the results computed with the two equations remain close to each other. The pressure distribution computed from Eq. (22), which employs the averaged channel height, is slightly higher than that from Eq. (21), which employs the effective channel height. However, the difference is no larger than 5% for almost the whole range of *λ* values considered. (The maximum deviation is found to be 5.14% in the case of *λ* = 10, which is pushing the limit of the applicability of the theory.) Having computed *P*(*Z*), $He(Z)$ and $H\xaf(Z)$ can be found from Eqs. (13) and (15), respectively. The largest deformed height is at the channel inlet (i.e., at *Z *=* *0), and we can expect the approximation of the effective channel height by the averaged one to be worst there. However, we determined that $max0\u2264\lambda \u226410|He(0)\u2212H\xaf(0)|/He(0)<5%$, showing good agreement.

Now that the validity of the approximate prediction of the flow rate–pressure drop relation given in Eq. (22) has been established, it is worthwhile to provide a formula for the fitting parameter *α* introduced by Gervais *et al.*^{22} Recall the averaged channel height from the latter model is

For a clearer comparison, we transform Eq. (15) into its dimensional form,

which we observe is a function of the Poisson's ratio, but no other material or geometric parameters related to the top wall, in this thick-wall limit ($t2/w2\u226b1$). [This observation will be contrasted with the result in Eq. (36) below.] Furthermore, most microchannels are made from materials such as polydimethylsiloxane (PDMS),^{42,43} which is often considered a nearly incompressible material, i.e., $\nu \u22480.5$. Then, $\alpha \u22480.4071$. A different solid mechanics model (and response) for the top wall would yield a different estimate of *α* (see Sec. III C 2), showing that *α* is *not* a universal number that can be determined by a single set of experiments (even if this approach works for some sets of geometries). Nevertheless, Eq. (29) provides a quantitative connection between the earlier scaling models^{22} for flow-induced deformation and the later detailed elasticity calculations.^{26}

It is also relevant to mention that the results in this subsection also yield insights into the quality of approximation of another approach to the flow-induced deformation problem. For example, following Skotheim and Mahadevan^{36,37} and Chakraborty and Chakraborty,^{38} Mukherjee *et al.*^{44} expressed the deformation at the fluid–solid interface of a thick-walled 2D duct as

where the layer thickness *H*_{1} and its “effective” Young's modulus *E _{m}* can be considered adjustable parameters.

^{45}In particular,

*H*

_{1}represents the distance over which the vertical displacement varies, vanishing at

*y*=

*H*

_{1}. Equation (30) is based on assuming no spanwise variation, reducing the flow and deformation problem to a 2D setting in the (

*y*,

*z*) plane, thus $h=h(z)$

*a forteriori*now (no averaging). The obvious question that arises is that what are suitable values of

*H*

_{1}and

*E*? As with Eq. (27), we simply compare Eq. (30) to Eq. (28) to obtain the answer. We conclude that

_{m}For example, if the 2D soft layer is taken to have the same elastic properties as the 3D one it approximates, *E _{m}* =

*E*, then Eq. (31) provides its suitable thickness

*H*

_{1}as a function of

*ν*and

*w*. Note that, separately, Essink

*et al.*

^{46}surveyed a number of such two-dimensional elastohydrodynamic problems, while Chandler and Vella

^{47}critically addressed the 2D models' validity in the near-incompressible limit as $\nu \u21921/2\u2212$.

#### 2. Duct with plate-like compliant top wall

Next, we analyze the case of a duct with a clamped thick-plate-like compliant top wall. As in Sec. III C 1, the slenderness of the duct still results in the decoupling of the top wall deformation at each streamwise cross section. However, the resulting shape of the deformed fluid–solid interface obtained by Shidhore and Christov^{29} is quite different from Eqs. (24) and (25). Specifically, now

where $B=E\xaft3/12$ is the plate's flexural rigidity,^{48} and *κ* is the “shear correction factor.”^{49} For consistency with the theory of elasticity, *κ* = 1 should be imposed,^{50} but we leave it in the equations for the sake of completeness. The plate model considers bending deformation, as well as shear deformation, of the top wall, and it is applicable for $t\u2272w$. If $t2/w2\u226a1$, the first term in the inner curly brace in Eq. (33) is negligible, meaning the shear deformation is not important in this case. The model then reduces to the one derived earlier by Christov *et al.*,^{23} which only accounted for plate bending.

By making Eq. (32) dimensionless, we obtain

Again, we have $F(X)\u2261F(x/w)=f(x)$, but $f$ is now given by Eq. (33). Then, Eq. (32) takes the same form as Eq. (16). Next, the calculation of the $Ii$ can be done explicitly for this case, yielding the functions of *t*/*w*, *κ* and *ν* summarized in Table II.

$I1$ | $130+(t/w)23\kappa (1\u2212\nu )$ |

$I2$ | $1630+(t/w)235\kappa (1\u2212\nu )+2(t/w)415[\kappa (1\u2212\nu )]2$ |

$I3$ | $112\u2009012+(t/w)2462\kappa (1\u2212\nu )+2(t/w)4105[\kappa (1\u2212\nu )]2+2(t/w)635[\kappa (1\u2212\nu )]3$ |

$I1$ | $130+(t/w)23\kappa (1\u2212\nu )$ |

$I2$ | $1630+(t/w)235\kappa (1\u2212\nu )+2(t/w)415[\kappa (1\u2212\nu )]2$ |

$I3$ | $112\u2009012+(t/w)2462\kappa (1\u2212\nu )+2(t/w)4105[\kappa (1\u2212\nu )]2+2(t/w)635[\kappa (1\u2212\nu )]3$ |

As in Sec. III C 1, we now substitute Eq. (34) into Eqs. (13) and (15), respectively, and compare the results. Figure 3 shows *P*(*Z*) for different *λ* and *Q *=* *1. The two formulations predict similar results. The error committed by replacing *H _{e}* with $H\xaf$ is $<8%$. However, even with smaller or larger

*t*/

*w*ratios, the pressure distributions computed with each

*H*expression do not differ much from each other. The maximum deviation is $<9%$. As in Sec. III C 1, we computed the absolute difference between using $He(0)$ and $H\xaf(0)$, and found that $max0\u2264\lambda \u226410|He(0)\u2212H\xaf(0)|/He(0)<5%$.

Finally, we can also compare the model from Eq. (22) (formulated with the averaged channel height) to Eq. (27) (the model derived by Gervais *et al.*^{22}) to obtain an explicit expression for the fitting parameter *α*. Again, for convenience, we write the dimensional form of the averaged channel height as

It follows, in this case, that

### D. A fiction factor for laminar flow in compliant ducts

Recently, it has been of interest to extend the textbook notion of a friction factor for various flows in microchannels. One idea is to take into account the shear-rate-dependent viscosity of non-Newtonian fluids.^{52} Even for Newtonian fluids, updates are being sought to better understand (the previously considered “settled”) wall roughness effects in both the laminar^{53} and turbulent^{54} portions of the Moody diagram (the visual representation of the friction factor^{55}) A friction factor is needed for microfluidic system design,^{56} much like its use for analyzing pipe networks.^{4} A frontier application is *microrheometry*,^{57,58} in which an experimentally computed friction factor in a rectangular microchannel is compared to a theoretical value, in order to characterize the viscosity of a fluid.^{59} An open problem in microrheometry^{60} concerns whether measurements made in PDMS microchannels are affected by the friction factor's implicit $\Delta p/E$ (or, in the present notation, *λ*) dependence. As the discussion above makes clear, the deformation of a compliant duct indeed changes the pressure drop characteristics. Thus, a salient application of our reduced-order flow and deformation model from Sec. III is to interrogate the dependence of the friction factor on the elasticity-related parameters and variables.

To this end, we start from the reduced model with the averaged channel height as the effective channel height, i.e., $he(z)=h\xaf(z)=h0[1+\eta p(z)]$. Note the compliance constant $\eta =\xi /Pc$, with $\xi =\lambda I1$ being the dimensionless spring stiffness parameter introduced in Sec. III A, is known from having solved a suitable solid mechanics problem. Then, from Eq. (22), we have

where we have substituted *Q *=* *1 and $Z=z/\u2113$. Equation (37) indicates that *ηp* cannot be varied independently because it is fully determined by *ξ*. In the following discussion, we work with dimensional variables for convenience.

For $Re\u0302\u21920$, the pressure difference across an axial length of a duct is balanced by the viscous drag on the wall. Denote the area of the cross section as $a(z)=wh\xaf(z)$, which takes into account the area change due to the deformation of the top wall. Then, the mean shear stress^{2} can be written as

Here, $cp=2(w+h\xaf)$ is the perimeter of the cross section, and $cp\u22482(w+h0)$ for $h\xaf\u226aw$. Additionally, $Dh0=4h0w/[2(w+h0)]$ is the *hydraulic diameter* of a rigid rectangular duct.^{2} In the last equality in Eq. (38), we further defined the hydraulic diameter of the soft duct as

where *ηp* captures the flow-induced deformation, meaning that *D _{h}* varies along the streamwise direction with

*p*.

Next, consider the Fanning friction factor defined^{2} as

where we have substituted Eq. (12) with $he=h\xaf$ into the last step above. Also note that we have introduced the averaged velocity as $v\xafz=q/(wh\xaf)$ and the hydraulic-diameter-based Reynolds number as

with $ReDh0=\rho qDh0/(\mu wh0)$ being the Reynolds number for the rigid rectangular duct.

Equation (40) has a form similar to the friction factor for a rigid rectangular duct. However, all three parameters, *D _{h}*, $h\xaf$, and $ReDh$, depend on

*z*due to FSI. To highlight this effect, we can re-write Eq. (40) as

The first term in Eq. (42) is *C _{f}* for a rigid rectangular duct, while the second term (in the parentheses) above captures the soft hydraulic effect. Furthermore, we can define the Poiseuille number as

We re-iterate that Eq. (43) is valid only for $h0\u226aw$ and observe that the prefactor $6(Dh0/h0)2=24$ for $h0/w\u21920$. Furthermore, while $Po=const.$ in a non-circular *rigid* duct,^{2} *Po* from Eq. (43) becomes a function of *z* due to FSI. Further, the increase in the soft hydraulic *Po* is clearly demonstrated by the second term in the last parenthesis on the right-hand side of Eq. (43), which is bounded between 1 (as $\eta p\u21920$) and 4 (as $\eta p\u2192\u221e$).

We highlight the novel dependence of *Po* on the compliance parameter *ξ*, beyond the usual geometric dependence on $(Dh0/h0)2$, by plotting *Po* vs $z/\u2113$ for given *ξ*, after eliminating *ηp* via Eq. (37). As predicted by Eq. (43), Fig. 4 shows that *Po* in a compliant duct is not a constant but rather a *decreasing* function along the streamwise direction [since *p*(*z*) is as well]. The shape is strongly influenced by the value of *ξ*, even if ultimately the correction factor due to compliance is bounded between 1 and 4.

## IV. SMALL BUT FINITE FLOW INERTIA: $Re\u0302=O(1)$

A feature of soft hydraulics problems is that the unidirectional flow solutions are derived under the lubrication approximation. As such, these solutions are *approximate* solutions and, thus, are *not* valid for arbitrary Reynolds number, unlike classical unidirectional flow solutions in ducts.^{2} Specifically, when the reduced Reynolds number, $Re\u0302$, is no longer vanishingly small, the inertial terms in Eq. (5) are no longer negligible. However, Eqs. (3) and (4) dictate that the pressure at each cross section is still uniform at the leading order (in *ϵ*); hence, we can still construct a 1D model relating the pressure *P*(*Z*) to the flow rate *Q*.

Toward this end, as before, we can either introduce $He(Z)$, based on enforcing a Poiseuille-like law from Eq. (20), or introduce the averaged channel height $H\xaf(Z)$ as an approximation to $He(Z)$ in the same relation. As shown in Sec. III for $Re\u0302\u21920$, using $H\xaf$ in place of *H _{e}* commits a controllable error, and both approaches lead to similar results (as long as the deformation gradient is small). Instead of treating both cases for $Re\u0302=O(1)$, we refer the reader to the work by Wang and Christov,

^{1}who implemented the calculation based on $He(Z)$. In this section, we construct a 1D model using $H\xaf(Z)$.

### A. Pressure distribution using an averaged deformed channel height

To accomplish this task, the von Kármán–Pohlhausen approximation (see Sec. 4–6.5 of White's book^{2}) is employed to enforce a shape of the streamwise velocity profile, $VZ2D$, so that the flow rate in the deformed fluid domain can be obtained.^{61–63} That is, we assume a dimensionless parabolic axial velocity profile $VZ2D$, which is related to the dimensionless volumetric flow rate *Q*, as

As discussed in Sec. III, a profile as in Eq. (44) is dictated by the Navier–Stokes equations for $Re\u0302\u21920$ (lubrication flow) and is generally valid for laminar flows.^{2} An implicit assumption for using the velocity profile in Eq. (44) for finite $Re\u0302$ is that flow inertia is weak: streamlines remain parallel and no recirculation occurs. Of course, this means that the theory developed herein is not valid in regimes in which transitional or turbulent flows occur. Indeed, the target application of our study is microfluidics, in which turbulent flows are not expected (or, generally possible),^{3,8} although laminar flow with $Re\u0302=O(1)$ can be achieved.^{64–66}

Observe that this expression, based on an equivalent 2D flow with $H\xaf$ as the effective channel height, does not depend (or require integration) over *X*. It should be noted that in the thin films literature^{67,68} inertial corrections to lubrication theory are also formulated, going to higher orders. Instead of assuming a parabolic velocity profile as in Eq. (44), a polynomial is used, and the coefficients are determined by incorporating the cross-sectional momentum equations, with their relevant boundary conditions, as well as the necessary corrections to the pressure. This approach is beyond the scope of the present work, however, as we consider wide channels ($\delta \u226a1$), and there is no cross-sectional flow components (*V _{X}* or

*V*) at the leading order

_{Y}^{23}in

*δ*and

*ϵ*.

Next, substituting $H\xaf$ from Eq. (18) into Eq. (45), we once again obtain a separable first-order ODE for *P*(*Z*). Imposing the outlet BC, $P(1)=0$, Eq. (45) integrates to

where $\xi =\lambda I1$ as above. As before, in the flow-controlled regime, *Q *=* *1, and $\Delta P$ is found implicitly from Eq. (46). Meanwhile, in the pressure-controlled regime, after enforcing $P(0)=1$, Eq. (46) becomes a quadratic in *Q*, and it has only one positive root,

This expression generalizes Eq. (23) and shows the dependence on $Re\u0302$ explicitly in the inertial flow.

Since Eq. (46) is a polynomial in *P*, we can invert it to find the pressure distribution in the duct. Importantly, we expect $dP/dZ<0$ strictly for all $X\u2208[0,1]$ because of the assumption of laminar flow. Since $P(1)=0$, *P*(*Z*)* *>* *0 for all $Z\u2208[0,1)$, which actually imposes an upper bound on the allowed values of $Re\u0302$ and *λ*. To prove this bound, the leading-order term of the left-hand side of Eq. (46) is calculated to be $(1\u22126Re\u0302\xi Q2/5)P$, as $Z\u21921\u2212$, while the right-hand side is positive. To ensure *P*(*Z*)* *>* *0 as $Z\u21921\u2212$, we must require that

### B. An extension and regularization via weak tension

At first glance, the restriction in Eq. (48) might be puzzling, but it actually ensures a continuous, and thus physical, pressure distribution and wall deformation at the leading order. Since the local deformed height is linearly proportional to the local pressure at the leading order, prominent local deformation can be expected for sufficiently inertial flows and/or sufficiently compliant ducts. In the case for which the restriction in Eq. (48) is violated, the local deformation can be so large that it cannot transition smoothly near to zero at the outlet [to satisfy the boundary condition $P(1)=0$, equivalently $H\xaf(1)=1$]. Thus, the solution from Eq. (46) breaks down for $Re\u0302$ values that violate the restriction in Eq. (48).

In deriving Eq. (46), we used Eq. (18), which is a leading-order solution (in *ϵ*) based on a plane strain configuration of the elastic wall's deformation field. This solution does not take into account the reaction forces imposed by connectors at the inlet and outlet of the duct. In this sense, we can think of the solid mechanics problem as being essentially a boundary layer problem. The Winkler-foundation-like mechanism [Eq. (18)] is dominant outside the boundary layers, while some other mechanism plays a role within thin (boundary) layers near *Z *=* *0, 1 to regularize the problem and account for the fact that the displacements in the vicinity of the inlet (or outlet) of the channel are usually restricted by external connections.

Since the bulging of the top wall unavoidably introduces stretching along *Z* in the solid, a simple extension of Eq. (18), which can circumvent the restriction from Eq. (48), can be achieved by introducing weak constant tension into the formulation.^{69} Note the tension has to be “weak” to ensure the dominance of the Winkler-foundation-like mechanism. Other regularization mechanisms are also possible. For example, in the setting of elastic structures on top of thin fluid films, Peng and Lister^{70} considered bending and gravity in addition to tension as regularization mechanisms. However, weak tension is arguably the simplest mechanism relevant to microchannels.

Then, we may write down a governing equation for the deformed channel height,^{71}

As motivated above, the dimensionless tension parameter $\theta 2\u226a1$. In this way, Eq. (18) is precisely the outer solution of Eq. (49) with $\theta 2=0$. To give a physical expression for $\theta 2$, we transform Eq. (49) back into dimensional form,

where *f _{t}* denotes the constant tension force per unit width ($N\u2009m\u22121$), and $K=Pc/(\xi h0)$ is the

*effective*stiffness of the interface ($Pa\u2009m\u22121$). Then, clearly, $\theta 2=\xi fth0/(Pc\u21132)$.

The tension *f _{t}* can arise from two physical scenarios. First,

*f*can arise due to stretching along

_{t}*z*, which is caused by the bulging of the fluid–solid interface. In this scenario,

*f*can be estimated by averaging the elongation of the fluid–solid interface along

_{t}*z*,

^{72}

Here, $t\u2605$ represents the *effective* thickness of the fluid–solid interface. If the wall is thin, we can take $t\u2605=t$. However, if the compliant wall is thick, the displacement decays away from the fluid–solid interface, as shown for the thick-walled case in Ref. 26. In this case, taking $t\u2605=t$ tends to overestimate the tension effect. Further considerations would be needed to estimate $t\u2605$ in this case, which is beyond the scope of the current work.

In the second scenario, *f _{t}* is provided by the pre-tension arising from external connectors. On the one hand, the pre-tension needs to be large enough so that the deformation induced stretch is negligible. On the other hand, the pre-tension needs to be small to ensure that $\theta 2\u226a1$. For the purposes of this paper, we focus on the effect of weak tension, which is consistent with our use of linear elasticity.

Next, taking $d/dZ$ of both sides of Eq. (49) and substituting into Eq. (45), we obtain a nonlinear ODE for $H\xaf(Z)$,

At the inlet and outlet, the top wall is restricted from moving, so the BCs for Eq. (52) are

where the BC in Eq. (54) is a restatement of the outlet BC $P(1)=0$ in terms of $H\xaf$ using Eqs. (49) and (53). Equations (52)–(54) constitute a nonlinear *two-point boundary-value problem.*^{73} As before, in the flow-controlled situation, *Q *=* *1 and Eq. (52) can be solved for $H\xaf(Z)$ subject to the BCs Eqs. (53) and (54). In the pressure-controlled situation, *Q* is found as an eigenvalue after imposing $P(0)=1$ on Eqs. (52)–(54).

Now, the restriction in Eq. (48) can be relaxed in the context of Eq. (52), in which the weak tension tends to restrain the wall deformation and, thus, regularizes the problem. Of course, the extent of regularization depends on the value of $\theta 2$. For example, if $\lambda =1.0$ and $I1=0.542\u2009754$ (for the thick-walled microchannel), then the upper bound of validity of the model is $Re\u0302\u22481.5$ for *θ* = 0. However, if $\theta 2=10\u22124$, Eq. (52) can be solved up to $Re\u0302\u22482.0$. If $\theta 2$ is further increased to $10\u22123$, then Eq. (52) can be solved up to $Re\u0302\u22483.0$. For such a large value of $Re\u0302$, one can interpret the breakdown of Eq. (52) as the breakdown of the lubrication theory and, potentially, as a sign that the “full” iNS equations need to be solved instead. Next, we illustrate these observations and explain how Eq. (52) was solved numerically.

### C. Illustrated examples

Depending on the top wall's geometry, *ξ* in Eqs. (45) and (52) will take different forms, such as the thick wall case and plate-like top wall case introduced in Secs. III C 1 and III C 2, respectively. To make our discussion general, instead of considering the two cases separately, we regard *ξ* and *θ* as characteristic system parameters and discuss the corresponding solutions of Eqs. (45) and (52) to illustrate the regularization introduced in Sec. IV B. Equation (45) can be solved in two steps. First, invert Eq. (46) to get *P*(*Z*). Second, substitute *P*(*Z*) into Eq. (18) to get $H\xaf$. As for Eq. (52), we use the solve_bvp routine from the SciPy stack^{74} to obtain a numerical solution of the nonlinear two-point boundary value problem. After obtaining $H\xaf(Z)$, Eq. (49) can be used to obtain *P*(*Z*).

First, we investigate the tension effect by varying $\theta 2$ in Eq. (45) while keeping *Q*, $Re\u0302$ and *ξ* fixed. As shown in Fig. 5(a), with $\theta 2\u226a1$, the solutions to Eqs. (45) and (52) do not differ much from each other along most of the domain $Z\u2208[0,1]$, as required by the dominance of the Winkler-foundation-like mechanism of deformation. Since Eq. (45) only satisfies $H\xaf(1)=1$, in principle, two boundary layers could be expected near *Z *=* *0 and *Z *=* *1, respectively, to fulfill the remaining boundary conditions from Eqs. (53) and (54). However, as we have discussed in Sec. IV A, with $\theta 2=0$, Eq. (46) indicates that *P* varies linearly with *Z* as $Z\u21921\u2212$. Since $H\xaf$ is linearly proportional to *P* at the leading order in *θ*, $H\xaf$ should be linear in *Z* as $Z\u21921\u2212$; hence, $d2H\xaf/dZ2\u21920$ as $Z\u21921\u2212$. In other words, the outer solution actually satisfies the boundary condition from Eq. (54). Therefore, there is no boundary layer located near *Z *=* *1. This fact can also be seen in Fig. 5(a), where the left boundary layer is prominent (becoming thicker as $\theta 2$ is increased), while the outer solution agrees well with the full numerical solution near *Z *=* *1 for all values of $\theta 2$ shown.

The effect of $\theta 2$ on *P*(*Z*) is shown in Fig. 5(b). The key takeaway from this plot is that, while Eq. (45) always predicts *P*(*Z*) to be a decreasing function of *Z*, a positive pressure gradient is observed near *Z *=* *0 in the numerical solution to Eq. (52) for all $\theta 2\u22600$ considered. The reason for this positive pressure gradient near the inlet is that, due to the restriction on the displacement at *Z *=* *0, the area of the cross section undergoes a sharp change near *Z *=* *0. Since the flow rate is fixed at steady state, the axial velocity has to quickly reduce near *Z *=* *0. The observed positive pressure gradient facilitates this deceleration of the flow.

Next, we address the effect of fluid inertia by varying $Re\u0302$. In this case, we fix *Q *=* *1, $\xi =0.5$, and $\theta 2=10\u22124$. As shown in Fig. 6(a), as $Re\u0302$ increases, larger deformation of the wall is observed. Also, the deformation gradient along *Z* is larger for higher $Re\u0302$ because the pressure displays a sharper decrease with the increase in $Re\u0302$, which can be clearly seen in Fig. 6(b). Notably, $dP/dZ>0$ is also observed for the three cases of $Re\u0302\u22600$, which can be explained as before. However, $dP/dZ$ remains negative in the case of $Re\u0302=0$. This is because, in this case of negligible fluid inertia, the deceleration of the flow near the inlet is not as large as the other cases; thus, the positive pressure gradient is not necessary. Finally, we mention that instead of solving Eq. (52) numerically, we are able to obtain a uniformly valid asymptotic solution for $H\xaf(Z)$ and *P*(*Z*) using the method of matched asymptotic expansions.^{75} In particular, for the special case of $Re\u0302=0$, we are able to obtain explicit formulas for both $H\xaf(Z)$ and *P*(*Z*). The details of this calculation are provided in Appendix A. The dashed curves in Figs. 6(a) and 6(b) demonstrate that these asymptotic solution [Eqs. (A13) and (A14) in Appendix A] agrees well with the numerical solution.

As a supplement to our discussion above, typical values of the dimensional and dimensionless variables of a microchannel with a thick top wall are summarized in Table III. Here, $t2/w2=16\u226b1$; thus, Eq. (26) and Table I from Sec. III C 1 are applicable. The steady responses of the system under different flow rates are calculated from Eq. (49) and tabulated in Table IV. With the increase in the flow rate, the pressure drop, the maximum pressure within the channel, and the maximum deformation of the interface increase. As we have discussed, when the flow inertia is small (smaller $Re\u0302$), the maximum pressure occurs at the inlet of the channel. However, if the flow inertia is prominent, there is a positive pressure gradient near the inlet and thus, the maximum pressure is “pushed” inwards, away from the inlet.

Name . | Variable . | Typical value . | Unit . |
---|---|---|---|

Channel's length | $\u2113$ | 1.0 | cm |

Channel's undeformed height | h_{0} | 25 | μm |

Channel's width | w | 500 | μm |

Top wall's thickness | t | 2.0 | mm |

Solid's Young's modulus | E | 1.5 | MPa |

Solid's Poisson's ratio | ν | 0.5 | ⋯ |

Fluid's dynamic viscosity | μ | $1.0\xd710\u22123$ | Pa s |

Fluid's density | ρ | $1.0\xd7103$ | kg m^{−3} |

Inlet flow rate | q | See Table IV | μl min^{−1} |

Tension force per unit width | f _{t} | 400 | N m^{−1} |

Characteristic velocity scale | $Vc=q/(wh0)$ | ⋯ | m s^{−1} |

Characteristic pressure scale | $Pc=\mu Vc/(\u03f5h0)$ | ⋯ | kPa |

Pressure drop | $\Delta p=p(z=0)$ | See Table IV | kPa |

Maximum pressure | $pmax=max0\u2264z\u2264\u2113p(z)$ | See Table IV | kPa |

Maximum channel's deformed height | $h\xafmax$ = $max0\u2264z\u2264\u2113h\xaf(z)$ | See Table IV | μm |

Channel's height-to-length aspect ratio | $\u03f5=h0/\u2113$ | 0.0025 | ⋯ |

Channel's height-to-width aspect ratio | $\delta =h0/w$ | 0.05 | ⋯ |

Reduced Reynolds number | $Re\u0302=\u03f5\rho q/(w\mu )$ | See Table IV | ⋯ |

Dimensionless spring stiffness | $\xi =\lambda I1$ ($\lambda =wPc/(h0E\xaf),\u2009I1=0.542754$) | See Table IV | ⋯ |

Tension coefficient | $\theta 2=fth0\xi /(Pc\u21132)=ftwI1(E\xaf\u21132)$ | $5.427\xd710\u22124$ | ⋯ |

Name . | Variable . | Typical value . | Unit . |
---|---|---|---|

Channel's length | $\u2113$ | 1.0 | cm |

Channel's undeformed height | h_{0} | 25 | μm |

Channel's width | w | 500 | μm |

Top wall's thickness | t | 2.0 | mm |

Solid's Young's modulus | E | 1.5 | MPa |

Solid's Poisson's ratio | ν | 0.5 | ⋯ |

Fluid's dynamic viscosity | μ | $1.0\xd710\u22123$ | Pa s |

Fluid's density | ρ | $1.0\xd7103$ | kg m^{−3} |

Inlet flow rate | q | See Table IV | μl min^{−1} |

Tension force per unit width | f _{t} | 400 | N m^{−1} |

Characteristic velocity scale | $Vc=q/(wh0)$ | ⋯ | m s^{−1} |

Characteristic pressure scale | $Pc=\mu Vc/(\u03f5h0)$ | ⋯ | kPa |

Pressure drop | $\Delta p=p(z=0)$ | See Table IV | kPa |

Maximum pressure | $pmax=max0\u2264z\u2264\u2113p(z)$ | See Table IV | kPa |

Maximum channel's deformed height | $h\xafmax$ = $max0\u2264z\u2264\u2113h\xaf(z)$ | See Table IV | μm |

Channel's height-to-length aspect ratio | $\u03f5=h0/\u2113$ | 0.0025 | ⋯ |

Channel's height-to-width aspect ratio | $\delta =h0/w$ | 0.05 | ⋯ |

Reduced Reynolds number | $Re\u0302=\u03f5\rho q/(w\mu )$ | See Table IV | ⋯ |

Dimensionless spring stiffness | $\xi =\lambda I1$ ($\lambda =wPc/(h0E\xaf),\u2009I1=0.542754$) | See Table IV | ⋯ |

Tension coefficient | $\theta 2=fth0\xi /(Pc\u21132)=ftwI1(E\xaf\u21132)$ | $5.427\xd710\u22124$ | ⋯ |

q (μl min^{−1})
. | $Re\u0302$ (–) . | ξ (–)
. | $\Delta p$ (kPa) . | $pmax$ (kPa) . | $h\xafmax$ (μm)
. |
---|---|---|---|---|---|

1500 | 0.125 | 0.1737 | 140.96 | 140.96 | 42.55 |

6000 | 0.5 | 0.6947 | 250.55 | 266.16 | 59.74 |

12 000 | 1.0 | 1.3895 | 258.27 | 366.90 | 73.61 |

q (μl min^{−1})
. | $Re\u0302$ (–) . | ξ (–)
. | $\Delta p$ (kPa) . | $pmax$ (kPa) . | $h\xafmax$ (μm)
. |
---|---|---|---|---|---|

1500 | 0.125 | 0.1737 | 140.96 | 140.96 | 42.55 |

6000 | 0.5 | 0.6947 | 250.55 | 266.16 | 59.74 |

12 000 | 1.0 | 1.3895 | 258.27 | 366.90 | 73.61 |

## V. CONCLUSION

In the spirit of Frank M. White's summary of unidirectional flows^{2} in non-circular ducts, we critically discussed weakly unidirectional flows (under a lubrication scaling) in compliant ducts of initially rectangular cross section for both the vanishing and the finite Reynolds number cases. In doing so, we contributed to the recently developed theory of *soft hydraulics*. Attention was paid to the hydraulic resistance of such conduits during steady viscous flow (i.e., the flow rate–pressure drop relations, which are now nonlinear). In particular, we derived 1D reduced models from 3D results on fluid–structure interaction. In doing so, we synthesized and unified a variety of previous models (some justified only by empirical considerations). This kind of reduction has been sought (and is of general interest^{62}) for practical design considerations of microfluidic systems,^{76–79} such as for calibrating optics-free non-contact measurement techniques.^{80}

For inertialess unidirectional flow in a compliant duct, the pressure varies nonlinearly along the streamwise direction due to the FSI between the viscous fluid flow and the compliant wall. Due to the slenderness and shallowness of the duct, we are able to relate the nonlinear pressure gradient $dp/dz$ to the flow rate *q* at steady state. By introducing the concept of an effective channel height, we recovered the form of the classical Poiseuille-like law and, at the same time, reduced the original 3D flow problem to an equivalent 2D one.

Although averaged deformed channel heights have been used in the literature, the validity of such models was not previously established. We found that the averaged channel height (14) can be a good approximation to the consistent effective height introduced in Eq. (11). This conclusion is important because the averaged-height models yield explicit flow rate–pressure drop relations and are easily compared to other geometries such as axisymmetric cases. Interestingly, we showed that the averaged channel height has a universal expression as $H\xaf(Z)=1+\xi P(Z)$, where $\xi =\lambda I1$, for *both* thick-walled and thinner, plate-like-walled top walls. Even though the formula for the dimensionless compliance coefficient *ξ* is different in the two cases, we have justified the observation (from the end of Sec. III A) that a wide and shallow microchannel's top wall behaves like a Winkler foundation,^{34,35} in which the averaged channel height is determined by the local pressure and a proportionality constant.

The reduction of the 3D FSI problem to a 1D model using the averaged height concept also allowed us to generalize the textbook concept of a friction factor^{2,4} to compliant ducts. We showed that the soft hydraulic system's Poiseuille number *Po* (product of the Fanning friction factor *C _{f}* and the Reynolds number) can be between 1 and 4 times larger than that for a rigid duct. Importantly, for the compliant duct, both

*C*and

_{f}*Po*depend on the streamwise coordinate due to the non-constant pressure gradient. This novel result extends the laminar portion of the Moody diagram, in which roughness is unimportant, via a new compliance parameter that is important in microfluidics.

Additionally, we showed how to incorporate weak but finite flow inertia in the previous $Re\u21920$ models. The finite-*Re* model breaks down beyond a certain value of the product of *Re* and a compliance parameter *λ*. Weak tension near the inlet and outlet of the reduced 1D model was introduced to regularize this breakdown and to obtain uniformly valid pressure distributions (in the sense of matched asymptotics).

The present results pave the way toward understanding more complex unsteady soft hydraulic phenomena. Specifically, with all this in mind, we would like to revisit and extend the linear stability results from Ref. 1 to the reduced-order models derived herein. This analysis could shed new insight on “ultrafast mixing” and multifold reduction of the critical Reynolds number recently observed in experiments on flow in compliant microchannels.^{81–83} Elastic walls (or coatings) have been shown to alter the turbulent boundary layer energy budget in channels;^{84,85} thus, the transition to turbulence in soft hydraulic systems^{83} is also expected to have nontrivial departures from the classical picture.

## DEDICATION

This paper is dedicated, with respect and admiration, to Professor Frank M. White on the occasion of his 88th anniversary.

## ACKNOWLEDGMENTS

I.C.C. thanks the Department of Mechanical Engineering at IIT Kharagpur for its hospitality during his visit there in December 2019, during which the idea for this work was conceived. Insightful discussions and input from J. Chakraborty and P. Karan on Winkler foundations and 2D models and from V. Anand and K. A. Flack on friction factors in pipes are kindly acknowledged.

I.C.C.'s visit to IIT Kharagpur and this research were enabled by the Scheme for Promotion of Academic and Research Collaboration (SPARC), a Government of India Initiative, under Project Code SPARC/2018–2019/P947/SL. Additionally, X.W. and I.C.C. were partially supported by the U.S. National Science Foundation under Grant No. CBET-1705637.

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this theoretical study. Python script files for generating the plots, which are based on the equations in the text, are openly available in the Purdue University Research Repository at https://dx.doi.org/10.4231/37PY-K896 and/or from the corresponding author upon reasonable request.

### APPENDIX A: MATCHED ASYMPTOTIC SOLUTION FOR THE 1D MODEL WITH WEAK TENSION

For $\theta 2\u226a1$, Eq. (52) subject to the BCs in Eqs. (53) and (54) represents a singular perturbation problem.^{75} The outer solution $H\xafo(Z)$, which satisfies $H\xafo(1)=1$, is found by setting $\theta 2=0$,

In the boundary layer near *Z *=* *0 (“left” boundary layer), we introduce the rescaled coordinate $Z\u0302=Z/\theta $. Denote the left inner solution as $H\xafl(Z\u0302)$. Then, in terms of these new variables, Eq. (52) is transformed into

At the leading order, the last term in Eq. (A2) is negligible, and we integrate once to obtain

Now, consider the behavior of Eq. (A3) in the phase plane $(H,F)$, where we have defined $H:=H\xafl$ and $F:=dH\xafl/dZ\u0302$; $Z\u0302$ parametrizes integral curves (i.e., solutions) in this plane. Equation (A3) becomes

Fixed points of the system (A4) and (A5) are such that the right-hand sides vanish. Although the expression for the fixed point $(H\u2605,F\u2605)$ with and $H\u2605>0$ and $F\u2605=0$ is lengthy, it can be found. The solution of Eq. (A3) as $Z\u0302\u2192\u221e$ and $dH\xafl/dZ\u0302\u21920$ should match the outer solution $H\xafo$ as $Z\u21920$. Therefore, $H\u2605$ must be chosen to be precisely $H\xafo(0)$, which is the positive real root of Eq. (A1) with *Z *=* *0. Consequently, without needing the explicit formula for $H\u2605$, we obtain

Now, the inner solution in the left boundary layer is the integral curve in the $(H,F)$ plane starting at $H=1$ and ending at $H=H\xafo(0)$. To construct this curve, multiply both sides of Eq. (A3) by $dH\xafl/dZ\u0302$ and obtain a first integral,

To ensure that $H\u2605=H\xafo(0)$ remains the desired fixed point of the ODE, the constant of integration must be

Then, Eq. (A7) can be rewritten as

Equation (A9) is separable, so its solution can be written as

where the positive root is taken because it is expected that $dH/dZ\u0302>0$ and thus, $H\xafo(0)>H$, in the boundary layer. Performing the integration in Eq. (A10) yields an implicit solution,

where $m=6Re\u0302\xi Q2/[5H\xafo(0)2]$. Observe that if the criterion in Eq. (48) is satisfied, then *m *<* *1 follows, which is required for the solution (A10) to exist. Therefore, the restriction (48) is needed to obtain a meaningful outer solution to Eq. (A1). In the case for which the criterion (48) is violated, this asymptotic analysis will break down, which suggests that tension is no longer a sufficiently small effect. In that case, we can solve Eq. (52) numerically.

Inverting Eq. (A11) to get an explicit expression for $H\xafl(Z\u0302)$ is nontrivial. However, for the special case of $Re\u0302=0$, Eq. (A10) immediately gives an explicit solution,

As for the right boundary, near *Z *=* *1, the ODE does not exhibit a boundary layer structure for $\theta 2\u21920$, as we discussed in Sec. IV A. This fact is also shown by Fig. 5(a), from which it is evident that the numerical solutions of the full ODE agree well with the leading-order outer solution (outside the left boundary layer), for any $\theta 2\u226a1$.

The composite solution is obtained after subtracting the common part between inner and outer solutions,

with $H\xafl$ and $H\xafo$ given (implicitly) by Eqs. (A11) and (A1), respectively. Equation (49) can be used to obtain the asymptotic solution for *P*. The leading-order terms are

where we have used Eq. (A3) to compute $d2H\xafl/dZ\u03022$.

For $Re\u0302=0$, using Eq. (A12), the composite solution can be explicitly written as

Substituting Eq. (A15) into Eq. (49) [or, setting $Re\u0302=0$ in Eq. (A14)], we obtain the matched asymptotic solution for the pressure distribution as well,

## References

^{4}the accepted definition that we shall employ: “[a] flow is classified as one-, two-, or three-dimensional depending on the number of space coordinates required to specify the velocity field” (p. 24).

*not*exact solutions of the incompressible Navier–Stokes equations.

^{86}

^{87,88}