Reduced models of unidirectional flows in compliant rectangular ducts at finite Reynolds number

As Frank M. White notes in Sec. 3–3.3 of his iconic Viscous Fluid Flow textbook,² “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³) 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,”² 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.⁴ 

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,¹² including fast progress in the last decade.¹³ While FSI topics such as aeroelasticity¹⁴ and moderate-Reynolds-number blood flow in large arteries⁶ are now quite classical, the mechanical interaction between slow viscous flows and compliant conduits¹⁵ 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²⁰ and, eventually, to one-dimensional (1D) models (that only involve axial, or streamwise, variations).²¹ 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.

Consider a soft-walled microchannel (initially a rectangular duct), which exhibits flow-induced deformation due to Newtonian fluid flow through it.²² Denote the channel's undeformed height, width, length, and top wall thickness by h₀, 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 shallow²³ 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 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)|−1/2<X<+1/2,0<Y<H(X,Z),0<Z<1}$⁠, in terms of the dimensionless variables. Here, $Rê$ is the modified Reynolds number, defined as $Rê=ϵ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.² 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 $Rê→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 varying²⁷ unidirectional flow solution,²⁸ 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 $Rê=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 $Rê→0$ 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

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

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,²³ 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,² 

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.²² and Hardy et al.,³³ 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¯$ (or $H¯$⁠) from Eq. (14) [or Eq. (15)] is not equal to h_e (or H_e) from Eq. (11) [or Eq. (13)]. Importantly, the averaging approach (introducing $h¯$ instead of h_e) 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 implies² 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

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,⁴⁰ and rollers.⁴¹ 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

and

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

Interestingly, observe that $H¯$ in Eq. (18) is simply the one-term Taylor-series approximation of H_e 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 Christov²⁶ and Shidhore and Christov²⁹) 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,

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¯$ 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 $Δ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.²² 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.,³³ 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).

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.²⁶ 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,

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.

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

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.²² Recall the averaged channel height from the latter model is

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

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

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 models²² for flow-induced deformation and the later detailed elasticity calculations.²⁶ 

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,³⁸ Mukherjee et al.⁴⁴ expressed the deformation at the fluid–solid interface of a thick-walled 2D duct as

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

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₁ as a function of ν and w. Note that, separately, Essink et al.⁴⁶ surveyed a number of such two-dimensional elastohydrodynamic problems, while Chandler and Vella⁴⁷ critically addressed the 2D models' validity in the near-incompressible limit as $ν→1/2−$⁠.

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

where $B=E¯t3/12$ is the plate's flexural rigidity,⁴⁸ and κ is the “shear correction factor.”⁴⁹ For consistency with the theory of elasticity, κ = 1 should be imposed,⁵⁰ 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.,²³ which only accounted for plate bending.

By making Eq. (32) dimensionless, we obtain

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.

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¯$ 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%$⁠.

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

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.⁵² Even for Newtonian fluids, updates are being sought to better understand (the previously considered “settled”) wall roughness effects in both the laminar⁵³ and turbulent⁵⁴ portions of the Moody diagram (the visual representation of the friction factor⁵⁵) A friction factor is needed for microfluidic system design,⁵⁶ much like its use for analyzing pipe networks.⁴ 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.⁵⁹ An open problem in microrheometry⁶⁰ 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

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 $Rê→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 stress² can be written as

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

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

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, D_h, $h¯$⁠, 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≪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,² 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.

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.² Specifically, when the reduced Reynolds number, $Rê$⁠, 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 $Rê→0$⁠, using $H¯$ 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 $Rê=O(1)$⁠, we refer the reader to the work by Wang and Christov,¹ 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 book²) 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 $Rê→0$ (lubrication flow) and is generally valid for laminar flows.² An implicit assumption for using the velocity profile in Eq. (44) for finite $Rê$ 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 $Rê=O(1)$ can be achieved.^64–66 

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

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 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 (⁠$δ≪1$⁠), and there is no cross-sectional flow components (V_X or V_Y) at the leading order²³ 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

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,

This expression generalizes Eq. (23) and shows the dependence on $Rê$ 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 $Rê$ and λ. To prove this bound, the leading-order term of the left-hand side of Eq. (46) is calculated to be $(1−6Rêξ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

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 $Rê$ 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.⁶⁹ 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⁷⁰ 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,⁷¹ 

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,

where f_t 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 f_t can arise from two physical scenarios. First, f_t can arise due to stretching along z, which is caused by the bulging of the fluid–solid interface. In this scenario, f_t can be estimated by averaging the elongation of the fluid–solid interface along z,⁷² 

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, 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 $θ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)$⁠,

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¯$ using Eqs. (49) and (53). Equations (52)–(54) constitute a nonlinear two-point boundary-value problem.⁷³ 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 $Rê≈1.5$ for θ = 0. However, if $θ2=10−4$⁠, Eq. (52) can be solved up to $Rê≈2.0$⁠. If $θ2$ is further increased to $10−3$⁠, then Eq. (52) can be solved up to $Rê≈3.0$⁠. For such a large value of $Rê$⁠, 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 stack⁷⁴ 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, $Rê$ 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.

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 $Rê$⁠. In this case, we fix Q = 1, $ξ=0.5$⁠, and $θ2=10−4$⁠. As shown in Fig. 6(a), as $Rê$ increases, larger deformation of the wall is observed. Also, the deformation gradient along Z is larger for higher $Rê$ because the pressure displays a sharper decrease with the increase in $Rê$⁠, which can be clearly seen in Fig. 6(b). Notably, $dP/dZ>0$ is also observed for the three cases of $Rê≠0$⁠, which can be explained as before. However, $dP/dZ$ remains negative in the case of $Rê=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.⁷⁵ In particular, for the special case of $Rê=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.