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.

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, Bruus3) 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, $−dp/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 biomechanics5–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 aeroelasticity14 and moderate-Reynolds-number blood flow in large arteries6 are now quite classical, the mechanical interaction between slow viscous flows and compliant conduits15 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) problems20 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.

FIG. 1.

Diagram of a long, shallow rectangular microchannel with a compliant top wall, labeled with the dimensional variables of the problem (denoted by lower case letters and symbols). The origin of the coordinate system is set at the centerline (x =0) of the rigid bottom wall of the channel (y =0). The deformed fluid–solid interface is defined as $y=h0+uy(x,z)$, where uy denotes the compliant top wall's y-displacement evaluated at y = h0. The Newtonian fluid flow, with a given volumetric flow rate q, is in the positive z-direction, as indicated by arrows, from the inlet at z =0 to the outlet at $z=ℓ$. Exemplar deformation profiles of the fluid–solid interface at different streamwise locations are shown by the red dashed curves, while the interface deformation along x =0 is represented by the red dash-dotted curve. Zero displacement conditions are enforced along z =0, $z=ℓ, x=−w/2$, and $x=w/2$.

FIG. 1.

Diagram of a long, shallow rectangular microchannel with a compliant top wall, labeled with the dimensional variables of the problem (denoted by lower case letters and symbols). The origin of the coordinate system is set at the centerline (x =0) of the rigid bottom wall of the channel (y =0). The deformed fluid–solid interface is defined as $y=h0+uy(x,z)$, where uy denotes the compliant top wall's y-displacement evaluated at y = h0. The Newtonian fluid flow, with a given volumetric flow rate q, is in the positive z-direction, as indicated by arrows, from the inlet at z =0 to the outlet at $z=ℓ$. Exemplar deformation profiles of the fluid–solid interface at different streamwise locations are shown by the red dashed curves, while the interface deformation along x =0 is represented by the red dash-dotted curve. Zero displacement conditions are enforced along z =0, $z=ℓ, x=−w/2$, and $x=w/2$.

Close modal

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 h0, w, $ℓ$, and t, respectively, as in Fig. 1. Further, introducing the aspect ratios, $δ=h0/w$ and $ϵ=h0/ℓ$, we say that the microchannel is long and shallow23 if $ϵ≪δ≪1$. 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)≪w≪ℓ$. 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 variables23,26 (denoted by capital letters),

$X=xw, Y=yh0, Z=zℓ,VX=δvxϵVc, VY=vyϵVc, VZ=vzVc, P=pPc.$
(1)

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

$∂VX∂X+∂VY∂Y+∂VZ∂Z=0,$
(2)
$−∂P∂X=0,$
(3)
$−∂P∂Y=0,$
(4)
$Rê(VX∂VZ∂X+VY∂VZ∂Y+VZ∂VZ∂Z)=−∂P∂Z+∂2VZ∂Y2.$
(5)

The fluid domain is defined as the deformed conduit: ${(X,Y,Z)|−1/2, in terms of the dimensionless variables. Here, $Rê$ is the modified Reynolds number, defined as $Rê=ϵRe=ϵρVch0/μ$, where ρ and μ are the fluid's density and dynamic viscosity, respectively. Equation (5) relates the characteristic pressure and velocity scales as $Pc=μℓVc/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 area2,3 via Eq. (5)] is not constant and also varies with Z.

The following discussion begins with the case of $Rê→0$ (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 varying27 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̂=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).

Neglecting the inertia of the flow by taking $Rê→0$ in Eq. (5), we find that the axial velocity VZ, subject to the no-slip boundary condition (BC) at the walls, has a parabolic variation along the height of the duct (Y-direction),

$VZ(X,Y,Z)=−12dPdZY[H(X,Z)−Y].$
(6)

At steady state, the flow rate is

$Q:=∫−1/2+1/2∫0H(X,Z)VZ(X,Y,Z) dY dX=const.,$
(7)

and thus the pressure gradient is found from Eqs. (6) and (7) to be

$−dPdZ=12Q∫−1/2+1/2H3(X,Z) dX.$
(8)

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=μqℓ/(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 $ΔP=P(Z=0)−P(Z=1)$ is controlled, then enforcing $P(0)=p(0)/Δp=1$ (i.e., taking $Pc=Δ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 $Δ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

$−dpdz=12μqwh03.$
(9)

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

$−dpdz=12μq∫−w/2+w/2h3(x,z) dx.$
(10)

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

$he(z):=[1w∫−w/2+w/2h3(x,z) dx]1/3.$
(11)

Then, Eq. (10) can be rewritten as

$−dpdz=12μqwhe3(z).$
(12)

Note that the corresponding dimensionless effective channel height is

$He(Z):=he(z)h0=[∫−1/2+1/2H3(X,Z) dX]1/3.$
(13)

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 he 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,

$h¯(z):=1w∫−w/2+w/2h(x,z) dx,$
(14)

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

$H¯(Z):=h¯(z)h0=∫−1/2+1/2H(X,Z) dX.$
(15)

It should be clear, however, that $h¯$ (or $H¯$) from Eq. (14) [or Eq. (15)] is not equal to he (or He) from Eq. (11) [or Eq. (13)]. Importantly, the averaging approach (introducing $h¯$ instead of he) leads to an inconsistency in the reduced model because if we replace $he3(z)$ with $h¯3(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.

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)≪w≪ℓ$] so that $uy≪w$ [recall that $h0≪w$ and $h(x,z)=h0+uy(x,z)$]. If the top wall is thick enough, with $w≲t≪ℓ$, 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≲w≪ℓ$, we require that $maxx,zuy≪t$ to make the linear elastic theory applicable.23,29,30 However, regardless of the wall thickness, as long as $t≪ℓ$, 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 implies2 that the shear stress $τyz≪p$. Since the tractions are continuous across the fluid–solid interface, we can thus infer that $σyz≪σyy$ where σyz and σyy 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 (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

$H(X,Z)=h0+uy(x,z)h0=1+λF(X)P(Z),$
(16)

where $λ:=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 $λ≪1/δ$ 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 $λ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¯(Z)$ as

$He(Z)=[1+3I1λP(Z)+3I2λ2P2(Z)+I3λ3P3(Z)]1/3$
(17)

and

$H¯(Z)=1+I1λP(Z).$
(18)

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

$Ii:=∫−1/2+1/2Fi(X) dX, i=1,2,….$
(19)

Interestingly, observe that $H¯$ in Eq. (18) is simply the one-term Taylor-series approximation of He from Eq. (17) in terms of $λ≪1$. However, our analysis does not require $λ≪1$; in fact, $λ=O(1)$ is possible. Linear elasticity only requires that $λ≪1/δ$ (as discussed by Wang and Christov26 and Shidhore and Christov29) Thus, we would like to determine if the approximation in going from Eq. (17) to Eq. (18) is a valid one.

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,

$−dPdZ=12QHe3(Z).$
(20)

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

$12Q(1−Z)=P(Z)[1+32I1λP(Z)+I2λ2P2(Z)+14I3λ3P3(Z)].$
(21)

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

$P(Z)=1I1λ{[48I1λQ(1−Z)+1]1/4−1}.$
(22)

As mentioned in Sec. III A, we may either consider a flow-controlled situation, in which Q =1 and $Δ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,

$Q=148×{4+6I1λ+4I2λ2+I3λ3, from (21),1I1λ[(I1λ+1)4−1], from (22).$
(23)

Equation (22) is essentially the same model derived by Gervais et al.22 However, in said work, $I1λ$ 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).

#### 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≳1.5$, the “thick” limit ($t2/w2≫1$) is achieved, and a simple analytical Fourier series solution can be written down for the deformed channel's top wall,

$h(x,z)=h0[1+wp(z)E¯h0f(x)],$
(24)
$f(x)=∑m=1∞2Ammπ sin [mπ(xw+12)],$
(25)

where we have defined $Am:=2mπ[1−(−1)m]$ and $E¯:=E/(1−ν2)$, with E being Young's modulus and ν the Poisson's ratio.

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

TABLE I.

Values of the coefficients ${Ii}i=13$ defined by Eq. (19) for the thick-walled microchannel.

 $I1$ 0.542 754 $I2$ 0.333 333 $I3$ 0.215 834
 $I1$ 0.542 754 $I2$ 0.333 333 $I3$ 0.215 834

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

$λ=wPch0E¯={μqℓh04E¯ (flow controlled),wΔph0E¯ (pressure controlled).$
(26)

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 $Δ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¯(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≤λ≤10|He(0)−H¯(0)|/He(0)<5%$, showing good agreement.

FIG. 2.

Thick top wall: axial pressure distribution P(Z) in a soft hydraulic conduit for Q =1 and different λ: (a) $λ=0.01$, (b) $λ=0.1$, (c) $λ=1.0$, and (d) λ = 10. The solid curve is computed from Eq. (21), in which the effective channel height from Eq. (13) is employed, while the dashed curve is computed from Eq. (22), in which the averaged channel height from Eq. (15) is employed. The shaded region represents $±5%$ of deviation from the solid curve, which is the baseline (or “truth”) for this model.

FIG. 2.

Thick top wall: axial pressure distribution P(Z) in a soft hydraulic conduit for Q =1 and different λ: (a) $λ=0.01$, (b) $λ=0.1$, (c) $λ=1.0$, and (d) λ = 10. The solid curve is computed from Eq. (21), in which the effective channel height from Eq. (13) is employed, while the dashed curve is computed from Eq. (22), in which the averaged channel height from Eq. (15) is employed. The shaded region represents $±5%$ of deviation from the solid curve, which is the baseline (or “truth”) for this model.

Close modal

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

$h¯(z)=h0[1+αwp(z)Eh0].$
(27)

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

$h¯(z)=h0[1+I1(1−ν2)wp(z)Eh0].$
(28)

Then, comparing Eqs. (27) and (28), it is readily recognized that

$α=I1(1−ν2)≈0.542 754(1−ν2),$
(29)

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≫1$). [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., $ν≈0.5$. Then, $α≈0.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 models22 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 Mahadevan36,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

$h(z)=h0[1+H1p(z)h0Em],$
(30)

where the layer thickness H1 and its “effective” Young's modulus Em can be considered adjustable parameters.45 In particular, H1 represents the distance over which the vertical displacement varies, vanishing at y = H1. 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 H1 and Em? As with Eq. (27), we simply compare Eq. (30) to Eq. (28) to obtain the answer. We conclude that

$H1Em=I1(1−ν2)wE≈0.542 754(1−ν2)wE.$
(31)

For example, if the 2D soft layer is taken to have the same elastic properties as the 3D one it approximates, Em = E, then Eq. (31) provides its suitable thickness H1 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 Vella47 critically addressed the 2D models' validity in the near-incompressible limit as $ν→1/2−$.

#### 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 Christov29 is quite different from Eqs. (24) and (25). Specifically, now

$h(x,z)=h0[1+w4p(z)24Bh0f(x)],$
(32)
$f(x)=[14−(xw)2]{2(t/w)2κ(1−ν)+[14−(xw)2]},$
(33)

where $B=E¯t3/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≲w$. If $t2/w2≪1$, 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

$λ=w4Pc24h0B={μqw3ℓ24h04B (flow controlled),w4Δp24h0B (pressure controlled).$
(34)

Again, we have $F(X)≡F(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.

TABLE II.

Functional forms of the coefficients ${Ii}i=13$ defined by Eq. (19) for the plate-like-walled microchannel.

 $I1$ $130+(t/w)23κ(1−ν)$ $I2$ $1630+(t/w)235κ(1−ν)+2(t/w)415[κ(1−ν)]2$ $I3$ $112 012+(t/w)2462κ(1−ν)+2(t/w)4105[κ(1−ν)]2+2(t/w)635[κ(1−ν)]3$
 $I1$ $130+(t/w)23κ(1−ν)$ $I2$ $1630+(t/w)235κ(1−ν)+2(t/w)415[κ(1−ν)]2$ $I3$ $112 012+(t/w)2462κ(1−ν)+2(t/w)4105[κ(1−ν)]2+2(t/w)635[κ(1−ν)]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 He with $H¯$ 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¯(0)$, and found that $max0≤λ≤10|He(0)−H¯(0)|/He(0)<5%$.

FIG. 3.

Plate-like top wall: Axial pressure distribution P(Z) in a soft hydraulic conduit for Q =1 and different λ: (a) $λ=0.01$, (b) $λ=0.1$, (c) $λ=1.0$, and (d) λ = 10. The solid curve is computed from Eq. (21), in which the effective channel height from Eq. (13) is employed, while the dashed curve is computed from Eq. (22), in which the averaged channel height from Eq. (15) is employed. The shaded region represents $±5%$ of deviation from the solid curve, which is the baseline (or “truth”) for this model. The top wall thickness-to-width ratio is $t/w=0.5$.

FIG. 3.

Plate-like top wall: Axial pressure distribution P(Z) in a soft hydraulic conduit for Q =1 and different λ: (a) $λ=0.01$, (b) $λ=0.1$, (c) $λ=1.0$, and (d) λ = 10. The solid curve is computed from Eq. (21), in which the effective channel height from Eq. (13) is employed, while the dashed curve is computed from Eq. (22), in which the averaged channel height from Eq. (15) is employed. The shaded region represents $±5%$ of deviation from the solid curve, which is the baseline (or “truth”) for this model. The top wall thickness-to-width ratio is $t/w=0.5$.

Close modal

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

$h¯(z)=h0[1+I1w4p(z)24Bh0]=h0[1+I1(1−ν22)(wt)3wp(z)Eh0].$
(35)

It follows, in this case, that

$α=I1(1−ν22)(wt)3=(1−ν260)[(wt)3+10κ(1−ν)(wt)].$
(36)

Observe that, unlike Eq. (29), α now depends upon w and t (with $w/t≳1$), in addition to ν. The dependence on t, which Eq. (36) now quantitatively predicts, has been observed in experimental studies.33,51

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 laminar53 and turbulent54 portions of the Moody diagram (the visual representation of the friction factor55) 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 microrheometry60 concerns whether measurements made in PDMS microchannels are affected by the friction factor's implicit $Δ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¯(z)=h0[1+ηp(z)]$. Note the compliance constant $η=ξ/Pc$, with $ξ=λ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

$ηp(z)=ξP(Z)=[48ξ(1−z/ℓ)+1]1/4−1,$
(37)

where we have substituted Q =1 and $Z=z/ℓ$. 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̂→0$, 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¯(z)$, which takes into account the area change due to the deformation of the top wall. Then, the mean shear stress2 can be written as

$τ¯w=−1cp(dpdza+pdadz)=−wh0cp(dpdz(1+ηp)+ηpdpdz)=Dh04(1+2ηp)(−dpdz)=Dh4(−dpdz).$
(38)

Here, $cp=2(w+h¯)$ is the perimeter of the cross section, and $cp≈2(w+h0)$ for $h¯≪w$. 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

$Dh:=Dh0(1+2ηp),$
(39)

where ηp captures the flow-induced deformation, meaning that Dh varies along the streamwise direction with p.

Next, consider the Fanning friction factor defined2 as

$Cf:=2τ¯wρv¯z2=12Dh2(−dpdz)(μρv¯zDh)(1μv¯z)=6(Dhh¯)21ReDh,$
(40)

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

$ReDh=ρv¯zDhμ=ReDh0(1+ηp1+ηp),$
(41)

with $ReDh0=ρqDh0/(μ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, Dh, $h¯$, and $ReDh$, depend on z due to FSI. To highlight this effect, we can re-write Eq. (40) as

$Cf=6(Dh0h0)21ReDh0︸rigid duct Cf(1+ηp1+ηp).$
(42)

The first term in Eq. (42) is Cf 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

$Po:=CfReDh=6(Dh0h0)2︸rigid duct Po(1+ηp1+ηp)2.$
(43)

We re-iterate that Eq. (43) is valid only for $h0≪w$ and observe that the prefactor $6(Dh0/h0)2=24$ for $h0/w→0$. 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 $ηp→0$) and 4 (as $ηp→∞$).

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/ℓ$ 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.

FIG. 4.

The variation of the reduced Poiseuille number $Po/(rigid duct Po)=[1+ηp/(1+ηp)]2$ from Eq. (43) along the flow wise direction, z, for different ξ.

FIG. 4.

The variation of the reduced Poiseuille number $Po/(rigid duct Po)=[1+ηp/(1+ηp)]2$ from Eq. (43) along the flow wise direction, z, for different ξ.

Close modal

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̂$, 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¯(Z)$ as an approximation to $He(Z)$ in the same relation. As shown in Sec. III for $Rê→0$, using $H¯$ in place of He 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̂=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¯(Z)$.

To accomplish this task, the von Kármán–Pohlhausen approximation (see Sec. 4–6.5 of White's book2) 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

$VZ2D(Y,Z)=6QY[H¯(Z)−Y]H¯3(Z).$
(44)

As discussed in Sec. III, a profile as in Eq. (44) is dictated by the Navier–Stokes equations for $Rê→0$ (lubrication flow) and is generally valid for laminar flows.2 An implicit assumption for using the velocity profile in Eq. (44) for finite $Rê$ 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̂=O(1)$ can be achieved.64–66

Substituting Eq. (44) into Eq. (5) and integrating over $Y∈[0,H¯(Z)]$, we obtain

$65RêddZ[Q2H¯(Z)]=−dPdZH¯(Z)−12QH¯2(Z).$
(45)

Observe that this expression, based on an equivalent 2D flow with $H¯$ as the effective channel height, does not depend (or require integration) over X. It should be noted that in the thin films literature67,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 ($δ≪1$), and there is no cross-sectional flow components (VX or VY) at the leading order23 in δ and ϵ.

Next, substituting $H¯$ 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

$P(Z)+32ξP2(Z)+ξ2P3(Z)+14ξ3P4(Z)−65RêξQ2P(Z)=12Q(1−Z),$
(46)

where $ξ=λI1$ as above. As before, in the flow-controlled regime, Q =1, and $Δ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,

$Q=25Rê2ξ2+524Rêξ(4+6ξ+4ξ2+ξ3)−5Rêξ.$
(47)

This expression generalizes Eq. (23) and shows the dependence on $Rê$ 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∈[0,1]$ because of the assumption of laminar flow. Since $P(1)=0$, P(Z)>0 for all $Z∈[0,1)$, which actually imposes an upper bound on the allowed values of $Rê$ and λ. To prove this bound, the leading-order term of the left-hand side of Eq. (46) is calculated to be $(1−6RêξQ2/5)P$, as $Z→1−$, while the right-hand side is positive. To ensure P(Z)>0 as $Z→1−$, we must require that

$Rêλ<56I1Q2.$
(48)

Note that $I1$ is set by the solution of the corresponding elasticity problem (recall Tables I and II).

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¯(1)=1$]. Thus, the solution from Eq. (46) breaks down for $Rê$ 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 Lister70 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

$−θ2d2H¯dZ2+H¯−1=ξP.$
(49)

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

$ftd2h¯dz2+K(h¯−h0)=p(z),$
(50)

where ft denotes the constant tension force per unit width ($N m−1$), and $K=Pc/(ξh0)$ is the effective stiffness of the interface ($Pa m−1$). Then, clearly, $θ2=ξfth0/(Pcℓ2)$.

The tension ft can arise from two physical scenarios. First, ft can arise due to stretching along z, which is caused by the bulging of the fluid–solid interface. In this scenario, ft can be estimated by averaging the elongation of the fluid–solid interface along z,72

$ft=Et*ℓ∫0ℓ12(dh¯dz)2 dz.$
(51)

Here, $t★$ represents the effective thickness of the fluid–solid interface. If the wall is thin, we can take $t★=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★=t$ tends to overestimate the tension effect. Further considerations would be needed to estimate $t★$ in this case, which is beyond the scope of the current work.

In the second scenario, ft 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 $θ2≪1$. 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¯(Z)$,

$35RêddZ(Q2H¯2)=1ξ(θ2d3H¯dZ3−dH¯dZ)−12QH¯3.$
(52)

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

$H¯(0)=H¯(1)=1,$
(53)
$d2H¯dZ2|Z=1=0,$
(54)

where the BC in Eq. (54) is a restatement of the outlet BC $P(1)=0$ in terms of $H¯$ 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¯(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 $θ2$. For example, if $λ=1.0$ and $I1=0.542 754$ (for the thick-walled microchannel), then the upper bound of validity of the model is $Rê≈1.5$ for θ = 0. However, if $θ2=10−4$, Eq. (52) can be solved up to $Rê≈2.0$. If $θ2$ is further increased to $10−3$, then Eq. (52) can be solved up to $Rê≈3.0$. For such a large value of $Rê$, 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.

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¯$. As for Eq. (52), we use the solve_bvp routine from the SciPy stack74 to obtain a numerical solution of the nonlinear two-point boundary value problem. After obtaining $H¯(Z)$, Eq. (49) can be used to obtain P(Z).

First, we investigate the tension effect by varying $θ2$ in Eq. (45) while keeping Q, $Rê$ and ξ fixed. As shown in Fig. 5(a), with $θ2≪1$, the solutions to Eqs. (45) and (52) do not differ much from each other along most of the domain $Z∈[0,1]$, as required by the dominance of the Winkler-foundation-like mechanism of deformation. Since Eq. (45) only satisfies $H¯(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 $θ2=0$, Eq. (46) indicates that P varies linearly with Z as $Z→1−$. Since $H¯$ is linearly proportional to P at the leading order in θ, $H¯$ should be linear in Z as $Z→1−$; hence, $d2H¯/dZ2→0$ as $Z→1−$. 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 $θ2$ is increased), while the outer solution agrees well with the full numerical solution near Z =1 for all values of $θ2$ shown.

FIG. 5.

(a) The deformed channel height $H¯(Z)$ for different values of the tension parameter $θ2$. The black curves represent the outer solution of Eq. (45), while the other (lighter) curves are obtained using the full (numerical) solution of the two-point boundary-value problem consisting of Eqs. (52)–(54). (b) The corresponding pressure distribution P(Z). The black curves are obtained by substituting the solution of Eq. (45) into Eq. (49), while the other (lighter) curves are similarly obtained from full (numerical) solution of Eqs. (52) and (53). In both panels, we fixed Q =1, $Rê=1.0$, and $ξ=0.5$.

FIG. 5.

(a) The deformed channel height $H¯(Z)$ for different values of the tension parameter $θ2$. The black curves represent the outer solution of Eq. (45), while the other (lighter) curves are obtained using the full (numerical) solution of the two-point boundary-value problem consisting of Eqs. (52)–(54). (b) The corresponding pressure distribution P(Z). The black curves are obtained by substituting the solution of Eq. (45) into Eq. (49), while the other (lighter) curves are similarly obtained from full (numerical) solution of Eqs. (52) and (53). In both panels, we fixed Q =1, $Rê=1.0$, and $ξ=0.5$.

Close modal

The effect of $θ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 $θ2≠0$ 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̂$. In this case, we fix Q =1, $ξ=0.5$, and $θ2=10−4$. As shown in Fig. 6(a), as $Rê$ increases, larger deformation of the wall is observed. Also, the deformation gradient along Z is larger for higher $Rê$ because the pressure displays a sharper decrease with the increase in $Rê$, which can be clearly seen in Fig. 6(b). Notably, $dP/dZ>0$ is also observed for the three cases of $Rê≠0$, which can be explained as before. However, $dP/dZ$ remains negative in the case of $Rê=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¯(Z)$ and P(Z) using the method of matched asymptotic expansions.75 In particular, for the special case of $Rê=0$, we are able to obtain explicit formulas for both $H¯(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.

FIG. 6.

(a) The deformed channel height $H¯(Z)$ for different values of the reduced Reynolds number $Rê$. The solid curves represent the numerical solution of Eq. (52), while the symbols represent the asymptotic solution [see Eqs. (A13) and (A15)]. (b) The corresponding pressure distributions P(Z). The solid curves are obtained by substituting the solution of Eq. (52) into Eq. (49), while the symbols are the asymptotic solution [see Eqs. (A14) and (A15)]. The agreement between the asymptotic and numerical solutions is so good that the curves mostly overlap. In both panels, we fixed Q =1, $ξ=0.5$, and $θ2=10−4$.

FIG. 6.

(a) The deformed channel height $H¯(Z)$ for different values of the reduced Reynolds number $Rê$. The solid curves represent the numerical solution of Eq. (52), while the symbols represent the asymptotic solution [see Eqs. (A13) and (A15)]. (b) The corresponding pressure distributions P(Z). The solid curves are obtained by substituting the solution of Eq. (52) into Eq. (49), while the symbols are the asymptotic solution [see Eqs. (A14) and (A15)]. The agreement between the asymptotic and numerical solutions is so good that the curves mostly overlap. In both panels, we fixed Q =1, $ξ=0.5$, and $θ2=10−4$.

Close modal

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≫1$; 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̂$), 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.

TABLE III.

Typical values of the dimensional and dimensionless parameters arising from Eq. (49).

NameVariableTypical valueUnit
Channel's length $ℓ$ 1.0 cm
Channel's undeformed height h0 25 μ
Channel's width w 500 μ
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×10−3$ Pa s
Fluid's density ρ $1.0×103$ kg m−3
Inlet flow rate q See Table IV  μl min−1
Tension force per unit width ft 400 N m−1
Characteristic velocity scale $Vc=q/(wh0)$ ⋯ m s−1
Characteristic pressure scale $Pc=μVc/(ϵh0)$ ⋯ kPa
Pressure drop $Δp=p(z=0)$ See Table IV  kPa
Maximum pressure $pmax=max0≤z≤ℓp(z)$ See Table IV  kPa
Maximum channel's deformed height $h¯max$ = $max0≤z≤ℓh¯(z)$ See Table IV  μ
Channel's height-to-length aspect ratio $ϵ=h0/ℓ$ 0.0025 ⋯
Channel's height-to-width aspect ratio $δ=h0/w$ 0.05 ⋯
Reduced Reynolds number $Rê=ϵρq/(wμ)$ See Table IV  ⋯
Dimensionless spring stiffness $ξ=λI1$ ($λ=wPc/(h0E¯), I1=0.542754$See Table IV  ⋯
Tension coefficient $θ2=fth0ξ/(Pcℓ2)=ftwI1(E¯ℓ2)$ $5.427×10−4$ ⋯
NameVariableTypical valueUnit
Channel's length $ℓ$ 1.0 cm
Channel's undeformed height h0 25 μ
Channel's width w 500 μ
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×10−3$ Pa s
Fluid's density ρ $1.0×103$ kg m−3
Inlet flow rate q See Table IV  μl min−1
Tension force per unit width ft 400 N m−1
Characteristic velocity scale $Vc=q/(wh0)$ ⋯ m s−1
Characteristic pressure scale $Pc=μVc/(ϵh0)$ ⋯ kPa
Pressure drop $Δp=p(z=0)$ See Table IV  kPa
Maximum pressure $pmax=max0≤z≤ℓp(z)$ See Table IV  kPa
Maximum channel's deformed height $h¯max$ = $max0≤z≤ℓh¯(z)$ See Table IV  μ
Channel's height-to-length aspect ratio $ϵ=h0/ℓ$ 0.0025 ⋯
Channel's height-to-width aspect ratio $δ=h0/w$ 0.05 ⋯
Reduced Reynolds number $Rê=ϵρq/(wμ)$ See Table IV  ⋯
Dimensionless spring stiffness $ξ=λI1$ ($λ=wPc/(h0E¯), I1=0.542754$See Table IV  ⋯
Tension coefficient $θ2=fth0ξ/(Pcℓ2)=ftwI1(E¯ℓ2)$ $5.427×10−4$ ⋯
TABLE IV.

Calculated steady-state responses of the microchannel system under different flow rates with the parameters specified in Table III.

q (μl min−1)$Rê$ (–)ξ (–)$Δp$ (kPa)$pmax$ (kPa)$h¯max$ (μ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̂$ (–)ξ (–)$Δp$ (kPa)$pmax$ (kPa)$h¯max$ (μ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

In the spirit of Frank M. White's summary of unidirectional flows2 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 interest62) 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¯(Z)=1+ξP(Z)$, where $ξ=λ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 factor2,4 to compliant ducts. We showed that the soft hydraulic system's Poiseuille number Po (product of the Fanning friction factor Cf and the Reynolds number) can be between 1 and 4 times larger than that for a rigid duct. Importantly, for the compliant duct, both Cf and 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→0$ 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 systems83 is also expected to have nontrivial departures from the classical picture.

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

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 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.

For $θ2≪1$, Eq. (52) subject to the BCs in Eqs. (53) and (54) represents a singular perturbation problem.75 The outer solution $H¯o(Z)$, which satisfies $H¯o(1)=1$, is found by setting $θ2=0$,

$1ξ[14(H¯o4−1)−65RêξQ2(H¯o−1)]=12Q(1−Z).$
(A1)

Substituting Eq. (18) into the above, we recover Eq. (46) as the outer solution for the pressure.

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

$35RêddẐ(Q2H¯l2)=1ξ(d3H¯ldẐ3−dH¯ldẐ)+θ12QH¯l3.$
(A2)

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

$35RêξQ2H¯l2=d2H¯ldẐ2−H¯l+C1.$
(A3)

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

$dHdẐ=F,$
(A4)
$dFdẐ=35RêξQ2H2+H−C1.$
(A5)

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

$C1=35RêξQ2H¯o(0)2+H¯o(0).$
(A6)

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¯o(0)$. To construct this curve, multiply both sides of Eq. (A3) by $dH¯l/dẐ$ and obtain a first integral,

$(dHdẐ)2=−65RêξQ2H+H2−2C1H+C2.$
(A7)

To ensure that $H★=H¯o(0)$ remains the desired fixed point of the ODE, the constant of integration must be

$C2=125RêξQ2H¯o(0)+H¯o(0)2.$
(A8)

Then, Eq. (A7) can be rewritten as

$(dHdẐ)2=[H−H¯o(0)]2{1−65RêξQ2HH¯o(0)2}.$
(A9)

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

$∫1H¯ldH[H¯o(0)−H]1−65RêξQ2HH¯o(0)2=Ẑ,$
(A10)

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

$−2[tanh−1(1−mH¯l)−tanh−1(1−m)]+21−mH¯o(0)[tanh−1(1−mH¯l1−mH¯o(0))−tanh−1(1−m1−mH¯o(0))]=Ẑ,$
(A11)

where $m=6RêξQ2/[5H¯o(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¯l(Ẑ)$ is nontrivial. However, for the special case of $Rê=0$, Eq. (A10) immediately gives an explicit solution,

$H¯l(Ẑ)=H¯o(0)+[1−H¯o(0)]e−Ẑ (Rê=0).$
(A12)

As for the right boundary, near Z =1, the ODE does not exhibit a boundary layer structure for $θ2→0$, 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 $θ2≪1$.

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

$H¯(Z)∼H¯a(Z)=H¯l(Z/θ)+H¯o(Z)−H¯o(0) (θ2≪1),$
(A13)

with $H¯l$ and $H¯o$ 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

$P(Z)∼Pa(Z)=1ξ(−d2H¯ldẐ2+H¯a−1)=1ξ{H¯o(Z)−1−35RêξQ2[1H¯l(Z/θ)2−1H¯o(0)2]}(θ2≪1),$
(A14)

where we have used Eq. (A3) to compute $d2H¯l/dẐ2$.

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

$H¯(Z)∼H¯a(Z)=[1−(1+48Qξ)1/4]e−Z/θ+[1+48Qξ(1−Z)]1/4 (θ2≪1, Rê=0).$
(A15)

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

$P(Z)∼Pa(Z)=1ξ{[1+48ξQ(1−Z)]1/4−1}(θ2≪1, Rê=0).$
(A16)
1.
X.
Wang
and
I. C.
Christov
, “
Soft hydraulics in channels with thick walls: The finite-Reynolds-number base state and its stability
,”
AIP Conf. Proc.
2302
,
020002
(
2020
).
2.
F. M.
White
,
Viscous Fluid Flow
, 3rd ed. (
McGraw-Hill Higher Education
,
New York
,
2006
).
3.
H.
Bruus
,
Theoretical microfluidics, Oxford Master Series in Condensed Matter Physics
(
Oxford University Press
,
Oxford, UK
,
2008
).
4.
P. J.
Pritchard
,
Fox & McDonald's Introduction to Fluid Mechanics
, 8th ed. (
John Wiley & Sons
,
Hoboken, NJ
,
2011
).
5.
S. I.
Rubinow
and
J. B.
Keller
, “
Flow of a viscous fluid through an elastic tube with applications to blood flow
,”
J. Theor. Biol.
35
,
299
313
(
1972
).
6.
T. J.
Pedley
,
The Fluid Mechanics of Large Blood Vessels
(
Cambridge University Press
,
Cambridge
,
1980
).
7.
J. B.
Grotberg
and
O. E.
Jensen
, “
Biofluid mechanics in flexible tubes
,”
Annu. Rev. Fluid Mech.
36
,
121
147
(
2004
).
8.
H. A.
Stone
,
A. D.
Stroock
, and
A.
Ajdari
, “
Engineering flows in small devices: Microfluidics toward a lab-on-a-chip
,”
Annu. Rev. Fluid Mech.
36
,
381
411
(
2004
).
9.
T. M.
Squires
and
S. R.
Quake
, “
Microfluidics: Fluid physics at the nanoliter scale
,”
Rev. Mod. Phys.
77
,
977
1026
(
2005
).
10.
G. M.
Whitesides
, “
The origins and the future of microfluidics
,”
Nature
442
,
368
373
(
2006
).
11.
Microfluidics and microscale transport processes, IIT Kharagpur Research Monograph Series
, edited by
S.
Chakraborty
(
CRC Press
,
Boca Raton, FL
,
2013
).
12.
M. P.
Païdoussis
,
Fluid-Structure Interactions: Slender Structures and Axial Flow
, Vol.
2
(
,
San Diego, CA
,
2016
).
13.
P.
Karan
,
J.
Chakraborty
, and
S.
Chakraborty
, “
Small-scale flow with deformable boundaries
,”
J. Indian Inst. Sci.
98
,
159
183
(
2018
).
14.
R. L.
Bisplinghoff
,
H.
Ashley
, and
R. L.
Halfman
,
Aeroelasticity
(
Dover Publications
,
Mineola, NY
,
1996
).
15.
D.
Chakraborty
,
J. R.
Prakash
,
J.
Friend
, and
L.
Yeo
, “
Fluid-structure interaction in deformable microchannels
,”
Phys. Fluids
24
,
102002
(
2012
).
16.
Fluid–Structure Interactions in Low-Reynolds-Number Flows
, edited by
C.
Duprat
and
H. A.
Stone
(
The Royal Society of Chemistry
,
Cambridge, UK
,
2016
).
17.
H.
Fallahi
,
J.
Zhang
,
H.-P.
Phan
, and
N.-T.
Nguyen
, “
Flexible microfluidics: Fundamentals, recent developments, and applications
,”
Micromachines
10
,
830
(
2019
).
18.
Y.
Matia
,
T.
Elimelech
, and
A. D.
Gat
, “
Leveraging internal viscous flow to extend the capabilities of beam-shaped soft robotic actuators
,”
Soft Rob.
4
,
126
134
(
2017
).
19.
P.
Polygerinos
,
N.
Correll
,
S. A.
Morin
,
B.
,
C. D.
Onal
,
K.
Petersen
,
M.
Cianchetti
,
M. T.
Tolley
, and
R. F.
Shepherd
, “
Soft robotics: Review of fluid-driven intrinsically soft devices; Manufacturing, sensing, control, and applications in human-robot interaction
,”
19
,
1700016
(
2017
).
20.
Note that these are 2D problems in an axial vertical (y, z) plane, using the axes notation in Fig. 1, not an (x, z) plane perpendicular to the flow direction as in White's 2D Dirichlet problems for unidirectional duct flows.2
21.
Since the issue of “dimensionality” of fluid flows and models has caused some confusion in the literature, here we restate. from Pritchard,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).
22.
T.
Gervais
,
J.
El-Ali
,
A.
Günther
, and
K. F.
Jensen
, “
Flow-induced deformation of shallow microfluidic channels
,”
Lab Chip
6
,
500
507
(
2006
).
23.
I. C.
Christov
,
V.
Cognet
,
T. C.
Shidhore
, and
H. A.
Stone
, “
Flow rate–pressure drop relation for deformable shallow microfluidic channels
,”
J. Fluid Mech.
841
,
267
286
(
2018
).
24.
Y.
Xia
and
G. M.
Whitesides
, “
Soft lithography
,”
Annu. Rev. Mater. Sci.
28
,
153
184
(
1998
).
25.
E.
Sollier
,
C.
Murray
,
P.
Maoddi
, and
D. D.
Carlo
, “
Rapid prototyping polymers for microfluidic devices and high pressure injections
,”
Lab Chip
11
,
3752
3765
(
2011
).
26.
X.
Wang
and
I. C.
Christov
, “
Theory of the flow-induced deformation of shallow compliant microchannels with thick walls
,”
Proc. R. Soc. A
475
,
20190513
(
2019
).
27.
M.
Van Dyke
, “
Slow variations in continuum mechanics
,”
25
,
1
45
(
1987
).
28.
Therefore, it is important to note that, unlike the case of unidirectional flows in rigid ducts, the solutions discussed herein are not exact solutions of the incompressible Navier–Stokes equations.86
29.
T. C.
Shidhore
and
I. C.
Christov
, “
Static response of deformable microchannels: A comparative modelling study
,”
J. Phys.: Condens. Matter
30
,
054002
(
2018
).
30.
V.
Anand
,
S. C.
Muchandimath
, and
I. C.
Christov
, “
Hydrodynamic bulge testing: Materials characterization without measuring deformation
,”
ASME J. Appl. Mech.
87
,
051012
(
2020
).
31.
S.
Rubin
,
A.
Tulchinsky
,
A. D.
Gat
, and
M.
Bercovici
, “
Elastic deformations driven by non-uniform lubrication flows
,”
J. Fluid Mech.
812
,
841
865
(
2017
).
32.
E.
Boyko
,
R.
Eshel
,
K.
Gommed
,
A. D.
Gat
, and
M.
Bercovici
, “
Elastohydrodynamics of a pre-stretched finite elastic sheet lubricated by a thin viscous film with application to microfluidic soft actuators
,”
J. Fluid Mech.
862
,
732
752
(
2019
).
33.
B. S.
Hardy
,
K.
Uechi
,
J.
Zhen
, and
H.
Pirouz Kavehpour
, “
The deformation of flexible PDMS microchannels under a pressure driven flow
,”
Lab Chip
9
,
935
938
(
2009
).
34.
E.
Winkler
,
Die Lehre von der Elastizität und Festigkeit mit besonderer Rücksicht auf ihre Anwendung in der Technik
(
Verlag von H. Dominicus
,
Prag
,
1867
).
35.
D. A.
Dillard
,
B.
Mukherjee
,
P.
Karnal
,
R. C.
Batra
, and
J.
Frechette
, “
A review of Winkler's foundation and its profound influence on adhesion and soft matter applications
,”
Soft Matter
14
,
3669
3683
(
2018
).
36.
J. M.
Skotheim
and
L.
, “
Soft lubrication
,”
Phys. Rev. Lett.
92
,
245509
(
2004
).
37.
J. M.
Skotheim
and
L.
, “
Soft lubrication: The elastohydrodynamics of nonconforming and conforming contacts
,”
Phys. Fluids
17
,
092101
(
2005
).
38.
J.
Chakraborty
and
S.
Chakraborty
, “
Influence of streaming potential on the elastic response of a compliant microfluidic substrate subjected to dynamic loading
,”
Phys. Fluids
22
,
122002
(
2010
).
39.
P.
Karan
,
J.
Chakraborty
, and
S.
Chakraborty
, “
Influence of non-hydrodynamic forces on the elastic response of an ultra-thin soft coating under fluid-mediated dynamic loading
,”
Phys. Fluids
32
,
022002
(
2020
).
40.
J.
Chakraborty
and
S.
Chakraborty
, “
Combined influence of streaming potential and substrate compliance on load capacity of a planar slider bearing
,”
Phys. Fluids
23
,
082004
(
2011
).
41.
X.
Yin
and
S.
Kumar
, “
Lubrication flow between a cavity and a flexible wall
,”
Phys. Fluids
17
,
063101
(
2005
).
42.
J. C.
McDonald
and
G. M.
Whitesides
, “
Poly(dimethylsiloxane) as a material for fabricating microfluidic devices
,”
Acc. Chem. Res.
35
,
491
499
(
2002
).
43.
J.
Friend
and
L.
Yeo
, “
Fabrication of microfluidic devices using polydimethylsiloxane
,”
Biomicrofluidics
4
,
026502
(
2010
).
44.
U.
Mukherjee
,
J.
Chakraborty
, and
S.
Chakraborty
, “
Relaxation characteristics of a compliant microfluidic channel under electroosmotic flow
,”
Soft Matter
9
,
1562
1569
(
2013
).
45.
Such models have been found useful in analyzing the global inflation or relaxation time scale of a microchannel, which is relevant to. the start-up problem and stop-flow lithography.87,88
46.
M. H.
Essink
,
A.
Pandey
,
S.
Karpitschka
,
C. H.
Venner
, and
J. H.
Snoeijer
, “
Regimes of soft lubrication
,”
J. Fluid Mech.
915
,
A49
(
2021
).
47.
T. G. J.
Chandler
and
D.
Vella
, “
Validity of Winkler's mattress model for thin elastomeric layers: Beyond Poisson's ratio
,”
Proc. R. Soc. A
476
,
20200551
(
2020
).
48.
S.
Timoshenko
and
S.
Woinowsky-Krieger
,
Theory of Plates and Shells
, 2nd ed. (
McGraw-Hill
,
New York
,
1959
).
49.
N.
Challamel
and
I.
Elishakoff
, “
A brief history of first-order shear-deformable beam and plate models
,”
Mech. Res. Commun.
102
,
103389
(
2019
).
50.
S.
Zhang
, “
On the accuracy of Reissner–Mindlin plate model for stress boundary conditions
,”
ESAIM: M2AN
40
,
269
294
(
2006
).
51.
M. K.
Raj
,
S.
DasGupta
, and
S.
Chakraborty
, “
Hydrodynamics in deformable microchannels
,”
Microfluid. Nanofluid.
21
,
70
(
2017
).
52.
Y. S.
Muzychka
and
J.
Edge
, “
Laminar non-Newtonian fluid flow in noncircular ducts and microchannels
,”
ASME J. Fluids Eng.
130
,
111201
(
2008
).
53.
Y.
Liu
,
J.
Li
, and
A. J.
Smits
, “
Roughness effects in laminar channel flow
,”
J. Fluid Mech.
876
,
1129
1145
(
2019
).
54.
K. A.
Flack
, “
Moving beyond Moody
,”
J. Fluid Mech.
842
,
1
4
(
2018
).
55.
F. L.
Moody
, “
Friction factors for pipe flow
,”
Trans. ASME
66
,
671
684
(
1944
).
56.
K. V.
Sharp
,
R. J.
,
J. G.
Santiago
, and
J. I.
Molho
, “
Liquid flows in microchannels
,” in
The MEMS Handbook
, edited by
M.
(
CRC Press
,
Boca Raton, FL
,
2001
), Chap. 6.
57.
C. J.
Pipe
and
G. H.
McKinley
, “
Microfluidic rheometry
,”
Mech. Res. Commun.
36
,
110
120
(
2009
).
58.
S.
Gupta
,
W. S.
Wang
, and
S. A.
Vanapalli
, “
Microfluidic viscometers for shear rheology of complex fluids and biofluids
,”
Biomicrofluidics
10
,
043402
(
2016
).
59.
X.
Yang
,
N. T.
,
Z.
Hayat
,
T.
Fu
,
S.
Jiang
,
C.
Zhu
, and
Y.
Ma
, “
Pressure drop of single phase flow in microchannels and its application in characterizing the apparent rheological property of fluids
,”
Microfluid. Nanofluid.
23
,
75
(
2019
).
60.
F.
Del Giudice
,
F.
Greco
,
P. A.
Netti
, and
P. L.
Maffettone
, “
Is microrheometry affected by channel deformation?
,”
Biomicrofluidics
10
,
043501
(
2016
).
61.
P. S.
Stewart
,
S. L.
Waters
, and
O. E.
Jensen
, “
Local and global instabilities of flow in a flexible-walled channel
,”
Eur. J. Mech. B
28
,
541
557
(
2009
).
62.
T. J.
Pedley
and
D.
Pihler-Puzović
, “
Flow and oscillations in collapsible tubes: Physiological applications and low-dimensional models
,”
Sādhāna
40
,
891
909
(
2015
).
63.
T. C.
Inamdar
,
X.
Wang
, and
I. C.
Christov
, “
Unsteady fluid-structure interactions in a soft-walled microchannel: A one-dimensional lubrication model for finite Reynolds number
,”
Phys. Rev. Fluids
5
,
064101
(
2020
).
64.
D.
Stoecklein
and
D.
Di Carlo
, “
Nonlinear microfluidics
,”
Anal. Chem.
91
,
296
314
(
2019
).
65.
D. D.
Carlo
,
D.
Irimia
,
R. G.
Tompkins
, and
M.
Toner
, “
Continuous inertial focusing, ordering, and separation of particles in microchannels
,”
Proc. Natl Acad. Sci. U. S. A.
104
,
18892
18897
(
2007
).
66.
E. J.
Lim
,
T. J.
Ober
,
J. F.
Edd
,
S. P.
Desai
,
D.
Neal
,
K. W.
Bong
,
P. S.
Doyle
,
G. H.
McKinley
, and
M.
Toner
, “
Inertio-elastic focusing of bioparticles in microchannels at high throughput
,”
Nat. Commun.
5
,
4120
(
2014
).
67.
C.
Ruyer-Quil
and
P.
Manneville
, “
Improved modeling of flows down inclined planes
,”
Eur. Phys. J. B
15
,
357
369
(
2000
).
68.
N. O.
Rojas
,
M.
Argentina
,
E.
Cerda
, and
E.
Tirapegui
, “
Inertial lubrication theory
,”
Phys. Rev. Lett.
104
,
187801
(
2010
).
69.
X. Y.
Luo
and
T. J.
Pedley
, “
A numerical simulation of unsteady flow in a two-dimensional collapsible channel
,”
J. Fluid Mech.
314
,
191
225
(
1996
).
70.
G. G.
Peng
and
J. R.
Lister
, “
Viscous flow under an elastic sheet
,”
J. Fluid Mech.
905
,
A30
(
2020
).
71.
Note that Eq. (49) would typically be written for the deformation, not the channel height. But, in our nondimensionalization, the deformation is simply $H¯−1$.
72.
I. J.
Hewitt
,
N. J.
Balmforth
, and
J. R.
De Bruyn
, “
Elastic-plated gravity currents
,”
Eur. J. Appl. Math.
26
,
1
31
(
2015
).
73.
H. B.
Keller
,
Numerical Solution of Two Point Boundary Value Problems, CBMS-NSF Regional Conference Series in Applied Mathematics
Vol.
24
(
SIAM
,
,
1976
).
74.
P.
Virtanen
,
R.
Gommers
,
T. E.
Oliphant
,
M.
Haberland
,
T.
Reddy
,
D.
Cournapeau
,
E.
Burovski
,
P.
Peterson
,
W.
Weckesser
,
J.
Bright
,
S. J.
van der Walt
,
M.
Brett
,
J.
Wilson
,
K. J.
Millman
,
N.
Mayorov
,
A. R. J.
Nelson
,
E.
Jones
,
R.
Kern
,
E.
Larson
,
C. J.
Carey
,
I.
Polat
,
Y.
Feng
,
E. W.
Moore
,
J.
VanderPlas
,
D.
Laxalde
,
J.
Perktold
,
R.
Cimrman
,
I.
Henriksen
,
E. A.
Quintero
,
C. R.
Harris
,
A. M.
Archibald
,
A. H.
Ribeiro
,
F.
Pedregosa
, and
P.
van Mulbregt
, “
SciPy 1.0: Fundamental algorithms for scientific computing in Python
,”
Nat. Methods
17
,
261
272
(
2020
).
75.
M. H.
Holmes
,
Introduction to Perturbation Methods
, Texts in Applied Mathematics (
Springer Science + Business Media
,
New York
,
2013
), Vol.
20
.
76.
A.
Mehboudi
and
J.
Yeom
, “
A one-dimensional model for compressible fluid flows through deformable microchannels
,”
Phys. Fluids
30
,
092003
(
2018
).
77.
A.
Mehboudi
and
J.
Yeom
, “
Experimental and theoretical investigation of a low-Reynolds-number flow through deformable shallow microchannels with ultra-low height-to-width aspect ratios
,”
Microfluid. Nanofluid.
23
,
66
(
2019
).
78.
H. A.
Stone
, “
Fundamentals of fluid dynamics with an introduction to the importance of interfaces
,” in
Soft Interfaces
, Lecture Notes of the Les Houches Summer School Vol.
98
, edited by
L.
Bocquet
,
D.
Quéré
,
T. A.
Witten
, and
L. F.
Cugliandolo
(
Oxford University Press
,
New York
,
2017
), pp.
3
79
.
79.
L.
Pekker
, “
A one-dimensional model of liquid laminar flows with large Reynolds numbers in tapered microchannels
,”
Phys. Fluids
33
,
042003
(
2021
).
80.
C.
Dhong
,
S. J.
Edmunds
,
J.
Ramírez
,
L. V.
Kayser
,
F.
Chen
,
J. V.
Jokerst
, and
D. J.
Lipomi
, “
Optics-free, non-contact measurements of fluids, bubbles, and particles in microchannels using metallic nano-islands on graphene
,”
Nano Lett.
18
,
5306
5311
(
2018
).
81.
M. K. S.
Verma
and
V.
Kumaran
, “
A multifold reduction in the transition Reynolds number, and ultra-fast mixing, in a micro-channel due to a dynamical instability induced by a soft wall
,”
J. Fluid Mech.
727
,
407
455
(
2013
).
82.
R.
Neelamegam
and
V.
Shankar
, “
Experimental study of the instability of laminar flow in a tube with deformable walls
,”
Phys. Fluids
27
,
024102
(
2015
).
83.
V.
Kumaran
and
P.
Bandaru
, “
Ultra-fast microfluidic mixing by soft-wall turbulence
,”
Chem. Eng. Sci.
149
,
156
168
(
2016
).
84.
M.
, “
Compliant coatings for drag reduction
,”
Prog. Aerosp. Sci.
38
,
77
99
(
2002
).
85.
M. E.
Rosti
and
L.
Brandt
, “
Low Reynolds number turbulent flows over elastic walls
,”
Phys. Fluids
32
,
083109
(
2020
).
86.
E.
Lauga
,
A. D.
Stroock
, and
H. A.
Stone
, “
Three-dimensional flows in slowly varying planar geometries
,”
Phys. Fluids
16
,
3051
3062
(
2004
).
87.
D.
Dendukuri
,
S. S.
Gu
,
D. C.
Pregibon
,
T. A.
Hatton
, and
P. S.
Doyle
, “
Stop-flow lithography in a microfluidic device
,”
Lab Chip
7
,
818
828
(
2007
).
88.
P.
Panda
,
K. P.
Yuet
,
D.
Dendukuri
,
T. A.
Hatton
, and
P. S.
Doyle
, “
Temporal response of an initially deflected PDMS channel
,”
New J. Phys.
11
,
115001
(
2009
).