Spontaneous oscillations and negative-conductance transitions in microfluidic networks

See allHide authors and affiliations

Science Advances  13 May 2020:
Vol. 6, no. 20, eaay6761
DOI: 10.1126/sciadv.aay6761


The tendency for flows in microfluidic systems to behave linearly poses challenges for designing integrated flow control schemes to carry out complex fluid processing tasks. This hindrance precipitated the use of numerous external control devices to manipulate flows, thereby thwarting the potential scalability and portability of lab-on-a-chip technology. Here, we devise a microfluidic network exhibiting nonlinear flow dynamics that enable new mechanisms for on-chip flow control. This network is shown to exhibit oscillatory output patterns, bistable flow states, hysteresis, signal amplification, and negative-conductance transitions, all without reliance on dedicated external control hardware, movable parts, flexible components, or oscillatory inputs. These dynamics arise from nonlinear fluid inertia effects in laminar flows that we amplify and harness through the design of the network geometry. These results, which are supported by theory and simulations, have the potential to inspire development of new built-in control capabilities, such as on-chip timing and synchronized flow patterns.


Microfluidic systems—networks of miniature flow channels capable of processing fluids—are now commonly used in applications ranging from chemical analysis (1) and flow cytometry (2) to computing (3) and point-of-care diagnostics (4). The value of microfluidic networks is manifest in their utility for manipulating fluid motion with precision. However, such manipulation is often controlled through the use of external hardware (46). For instance, microscopic valves generally need to be actuated by macroscopic, computer-operated pumps (7), which has impeded development of portable microfluidic systems (4, 6). The need for active control stems from the low Reynolds numbers typical of microfluidic flows, whereby fluid inertia forces are small relative to viscous dissipation, causing flow rate changes to be linearly related to pressure changes (8). Thus, it remains challenging to design integrated control mechanisms that are capable of inducing responsive flow dynamics, such as oscillations, switching, and amplification, without relying on nonlinear input signals or moveable parts.

Nonetheless, substantial progress has been made in the development of built-in microfluidic controls. State-of-the-art approaches for incorporating passive valves for flow-rate regulation generally take advantage of flexible membranes and surfaces to generate nonlinear fluid-structure interactions (911). Complex flow patterns and operations have been implemented in such networks, but flexible components can hinder integration, yield to high driving pressures, and may require polymer materials that are not chemically compatible with the working fluid (4, 6, 12). On the other hand, recent appreciation has emerged for the impact and utility of fluid inertia effects on manipulating local flow dynamics in microfluidics (13, 14). It has been shown that even for moderate Reynolds numbers, the formation of vortices and secondary flows can be exploited for particle segregation (1518), mixing fluids (19, 20), and diverting flow streams (21, 22).

Here, we present a microfluidic network construction that demonstrates new dynamics resulting from fluid inertia, which can serve as novel flow control mechanisms and facilitate the design of integrated microfluidic systems. Our network exhibits (i) spontaneous emergence of persistent flow-rate oscillations for fixed driving pressures; (ii) hysteretic flow behavior in which more than one set of stable flow rates exist for the same driving pressures; and (iii) negative-conductance transitions, whereby an increase (decrease) in the driving pressure leads to a discontinuous decrease (increase) in the flow rate. These behaviors are interesting in their own right and are analogous to behaviors formerly sought through different approaches. Oscillations have been implemented in microfluidic networks by using flexible components (23, 24) and used as a timing mechanism (11, 25). Hysteresis has been explored through the implementation of hysteretic valves and, along with oscillatory driving, found applications in establishing microfluidic logic systems (10, 26). Nonmonotonic pressure-flow relationships analogous to negative-conductance transitions have been previously sought using flexible diaphragm valves (27, 28) and used for signal amplification and flow switching. Our network does not include flexible components, nor does it rely on oscillatory inputs. Instead, the behaviors in (i) to (iii) arise by structuring the network so that dynamic vortices are generated in the flow and nonlinear fluid inertia effects are amplified. The results presented in this work are derived from simulations of the Navier-Stokes equations and an analytical dynamical model developed to capture the diverse flow properties of the network.

Microfluidic network description and simulation results

A circuit schematic of our microfluidic network is shown in Fig. 1. The network consists of five channel segments that are constructed into two parallel paths connected by a transversal path. Typically, the steady-state relation between the flow rate Q through a microfluidic channel and the pressure loss ΔP along the channel takes the form ΔP = RQ, where R is the (absolute) fluidic resistance of the channel. When R is constant, this relation is analogous to Ohm’s law for electronic resistors (29, 30). Therefore, we represent the three (straight) channels in the network that exhibit constant fluidic resistance as linear resistors in the schematic (Fig. 1A). The two remaining channels include either a chicane of blade-like barriers (Fig. 1B) or an array of six cylindrical obstacles (Fig. 1C) that induce nonlinear pressure-flow relations and are represented as nonlinear resistors. As we show below, the obstacle-laden channel serves to amplify inertial effects and the chicane channel gives rise to oscillations. The lengths of the channels vary (Fig. 1D) but all share a common width w of 500 μm. The cylindrical obstacles have a radius of w/5, and the two barriers, which extend to the center of the chicane channel, are of thickness w/10. No-slip boundary conditions are assumed at all surfaces, and we consider the static pressure at the outlets of the system Pout to be held at a fixed common value, taken to be zero. At the inlets, we control either the pressures (P1in and P2in) or the flow rates (Q1 and Q2).

Fig. 1 Microfluidic network structure.

(A) Circuit schematic of the network, where the labels denote pressures (Pi), channel resistances (Ri), and flow rates (Qi). The inlet and outlet pressures are identified by the superscripts “in” and “out”, respectively, and the positive flow directions are indicated by arrows. Two channels exhibit variable (flow-dependent) resistance due to the presence of obstacles. (B and C) Geometric structure of the chicane (B) and obstacle-laden (C) channels. The blue curves mark example streamlines and specify flow direction. The closed streamlines in (B) represent vortices that form near the barriers and r marks the linear size of the left vortex. (D) Network topology of the circuit in (A), where the length of each channel segment is labeled by Li.

We present the outstanding properties of this microfluidic network through fluid dynamics simulations of incompressible flow in two dimensions. We consider a water-like working fluid with density ρ = 1000 kg/m3 and dynamic viscosity μ = 10−3 Pa·s. In microfluidics, pressure-driven flow is used across a variety of applications (5), whereby the system inlets are connected to a pressurized fluid reservoir, the outlets are open to atmosphere (or a lower-pressure reservoir), and flow is driven by the resulting pressure gradient. Here, we investigate the case in which a common static pressure is applied at the inlets, that is, P1in=P2in=Pin, which corresponds to the physical scenario in which the inlets are connected to a high-pressure reservoir through intermediate passive pressure regulators.

In Fig. 2A, we show simulation results of the total flow rate QT = Q1 + Q2 through the network in Fig. 1 over a range of driving pressures, Pin, from which we observe two notable properties. First, for Pin within two disjoint ranges, two stable solutions for the total flow rate exist. Second, a subset of solutions are unsteady and exhibit oscillating flow rates (Supplementary Materials, fig. S3, and movie S1), despite Pin being fixed. In particular, we find that at a critical value of Pin, solutions along the high-flow branch (red symbols in Fig. 2A) become small-amplitude limit cycles. The corresponding amplitudes and periods grow with Pin (the frequency of the oscillations decreases from 20 to 4 Hz; see the Supplementary Materials and fig. S2). At a higher critical Pin, the limit cycle collides with the unstable branch, thereby destabilizing the high-flow solution branch. An important property of the oscillating solutions is that the proportions of the flow rates through different channel segments also become time dependent (Fig. 2B). Bistability and spontaneous oscillations have been previously studied in fixed-structure microfluidic networks when feedback loops are incorporated (31) or when multiple working fluids with different viscosities are used (32). However, neither of these mechanisms are required in our system.

Fig. 2 Bistability and spontaneous oscillations.

(A) Bifurcation diagram of total flow rate as a function of the inlet pressure Pin, generated from direct simulations of the network in Fig. 1 for P1in=P2in=Pin. There exist stable high-flow (red) and low-flow (blue) solution branches, separated by an unstable intermediate-flow branch (black). Oscillating solutions arise spontaneously along the high-flow branch, where the oscillation amplitude is indicated by the shaded region. The Reynolds number for flows through the chicane channel and the obstacle-laden channel are in the range of 14 to 90 and 80 to 155, respectively. The solutions for Pin = 180 Pa, marked with a, b, and c, will be used as references in comparing with other figures. The unstable solutions are determined through flow-controlled simulations (see Supplementary Materials and fig. S1). (B) Time series of the proportion of flow exiting the obstacle-laden channel that passes through the chicane channel (Q3/Q4) for two driving pressures that yield oscillatory flows, showing that frequency decreases and amplitude increases as the driving pressure is increased.

Another outstanding property that arises from the bistability in our system is the possibility of negative-conductance transitions and other sudden transitions in QT that result from small changes in Pin. We characterize these transitions, which occur at the boundaries of the bistable regions (Fig. 2A), by defining (local) fluidic conductance and resistance as C = δQTPin and its reciprocal, respectively. Here, δ indicates a finite change and Pin is the controlled variable. Therefore, negative-conductance and negative-resistance transitions occur when an increase (decrease) in Pin leads to a decrease (increase) in QT. Notably, as shown in Fig. 3A, our system exhibits transition points at which CPin) diverges in the limit of small δPin: two points at which CPin → 0) = +∞, corresponding to positive-conductance transitions, and two points at which CPin → 0) = −∞, corresponding to negative-conductance transitions. Figure 3B shows that related transitions emerge when the flow rate QT, rather than the pressure Pin, is taken as the control variable. In this case, a change in QT can lead to transitions in which Pin changes by a finite amount. In particular, the latter includes signal amplification transitions, which are remarkable transitions in which an infinitesimal increase (decrease) in QT leads to a finite decrease (increase) in Pin. Both the pressure- and flow-driven transitions reported here are intimately related to the emergence of hysteresis in the system, which is another consequence of bistability that has potential applications in the development of systems with built-in memory.

Fig. 3 Hysteresis and flow state transitions.

(A) Hysteresis loop and resulting negative-conductance transitions for the network in Fig. 1 when quasistatically increasing (red) or decreasing (blue) the inlet driving pressure. (B) Counterpart of (A) and resulting signal amplification transitions when quasistatically varying the total flow rate. For the latter, Q1 and Q2 are controlled to maintain equal pressures at the inlets.

The solutions belonging to the different branches in Fig. 2A can be further distinguished by the flow rate and internal flow structure through specific channels. It is particularly insightful to examine the streamlines around the complex geometry in the chicane channel and the associated flow rate Q3. In Fig. 4A, we show the streamlines corresponding to the three labeled states in Fig. 2A. A number of steady vortices are observed in the flow around the barriers. The sizes of the vortices are correlated and we designate r to be the size of one of them, as labeled in Figs. 1B and 4A. We use a one-dimensional measure for r, taken to be the distance from the barrier to the vortex reattachment point along the channel wall. In Fig. 4 (B and C), we show that both Q3 and r differ markedly for solutions belonging to the three branches in Fig. 2A and that oscillations simultaneously emerge in these variables (Supplementary Materials, fig. S3, and movie S1). Notably, solutions along the (high-) low-flow branch in Fig. 2A correspond to (large) small values of Q3 and r.

Fig. 4 Flow structure in the chicane channel.

(A) Streamlines corresponding to the labeled solutions in Fig. 2A show variations in the vortices around the blade barriers. The size of one of the vortices is denoted by r. (B and C) Bifurcation diagrams for Q3 (B) and r (C), corresponding to all simulation results presented in Fig. 2A.

We determine the relationship between r and Q3 by performing simulations in which the flow rates at both inlets (Q1 and Q2) are controlled. From these simulations, we compute r, Q3, and the pressure loss along the chicane channel ΔP34, where the latter corresponds approximately to P3P4 (Fig. 1). In Fig. 5, we show relations between these quantities for sets of simulations in which Q1 is fixed while Q2 is varied. We observe nonlinear relations between r and Q3 (Fig. 5A), between Q3 and the pressure loss along the chicane channel (Fig. 5B), and between r and the fluidic resistance of the chicane channel (Fig. 5C). These nonlinear relations suggest a coupling between the pressure-flow relation of the chicane channel and the vortex size. We also note that discontinuities arise in the pressure-flow relation for the chicane channel (Fig. 5B) that result from abrupt changes in the vortex size as Q2 is varied (Fig. 5A). These discontinuities show the emergence of regions where the pressure-flow relation is negatively sloped, which correspond to regions of negative differential resistance.

Fig. 5 Vortex-flow rate interaction.

(A to C) Navier-Stokes simulation results of the network in Fig. 1 for fixed values of Q1 as Q2 is increased, from which we determine the relation between Q3 and r (A), the pressure-flow relation for the chicane channel (B), and the dependence of the chicane channel resistance on r (C). Transitions are evident at the points of discontinuity in (B), which can be associated with the points of discontinuity in Fig. 3, albeit for different control and independent variables. The pressures P3 and P4 are approximated from simulations by averaging pressure values sampled across the channel width near the chicane channel junctions. The chicane channel resistance is defined as (P3P4)/Q3 and is nondimensionalized by dividing it by μ/w2.

Analytical dynamical model

We now construct an analytical model of the system in Fig. 1 that characterizes our simulation results. For unidirectional laminar flow through a straight channel, the average flow rate of an incompressible fluid can be approximated from the Navier-Stokes equations aslQ·=ΔPRQ(1)where the dot implies a time derivative and l may be referred to as the fluidic inductance (33). For flow through a two-dimensional channel of length L, where the characteristic time scale of the flow is larger than the viscous time scale, the fluidic resistance and inductance can be approximated as R = 12μL/w3 and l = ρL/w. More generally, when the time scale of the flow exceeds the viscous time, memory effects in R and l become substantial. Under steady flow conditions, Eq. 1 reduces to ΔP = RQ.

One of the assumptions in the derivation of Eq. 1 is that all streamlines of the channel flow are straight, which causes the nonlinear inertial terms in the Navier-Stokes equations to vanish. Streamlines in the chicane channel violate this assumption (Fig. 4A) and nonlinear effects are therefore expected to be present. Indeed, we observe an approximately quadratic relation between the chicane channel resistance, R3, and the vortex size, r, for 60 < r < 400 μm (Fig. 5C). To construct an approximate dynamical equation for Q3, we use the form of Eq. 1 with the constant resistance replaced by a function of r. Specifically, we take R3(r) = 12μ(Lb + γ(rrb)2)/w3, where Lb serves as a base component of the resistance, γ is a constant coefficient of the variable component that depends on the vortex size, and rb is the vortex size that minimizes the resistance (from Fig. 5C, rb ≈ 150 μm). With this added dependence on r, we must also account for the dynamics of the vortex size. The steady-state relation between Q3 and r found through flow-controlled simulations (Fig. 5A) can be well fit by a cubic equation of the form Q3Q3*=η(rr*)3ξ(rr*), where η and ξ are positive parameters and Q3* and r* are constants that shift the cubic relation from the origin. For simplicity, we consider the growth rate of r to be proportional to the deviation from this equilibrium relation. Therefore, the dynamical equations that characterize the chicane channel take the formr·=ε(Q3Q3*η(rr*)3+ξ(rr*))(2)Q·3=wΔP34ρL312νw2L3(Lb+γ(rrb)2)Q3(3)where ε is a positive constant. For suitable parameters, we find that these equations capture the most salient properties in Fig. 4 (B and C). We show in Fig. 6 that for different ΔP34, Eqs. 2 and 3 can exhibit bistability and stable limit cycle solutions. We note that the additional dependence of the relations presented in Fig. 5 on Q1 can be accounted for by allowing η and Lb to be functions of Q1 (Supplementary Materials).

Fig. 6 Analytical dynamical model of flow through the chicane channel.

(A to D) Phase space plots showing example trajectories (red and green curves) and streamlines (gray curves), generated from Eqs. 2 and 3 for the flow rate and vortex dynamics at different values of ΔP34. Fixed-point solutions to the equations exist at the intersections of the r and Q3 nullclines (i.e., the curves in phase space for which r·=Q·3=0). The solution set may consist of one steady solution (A), three steady solutions (two stable and one unstable) (B), two steady solutions (one stable and one unstable) and a stable limit cycle (C), or a single stable limit cycle (D), depending on the value of ∆P34. (E) Bifurcation diagram of Q3 produced for Eqs. 2 and 3. The parameters used here are Q3*=25μl/s per mm depth, rb = 146 μm, r* = 250 μm, Lb = 0.146 cm, γ = 0.264 μm−1, and ε = 4166 m−1.

A second nonlinear element of the network in Fig. 1 is the obstacle-laden channel. As the flow rate through this channel segment increases, stationary eddies form in the wake of the obstacles for moderate Reynolds numbers. The presence of many of these obstacles in close proximity generates large velocity gradients in the surrounding flow, which amplifies energy dissipation and results in an overall nonlinear pressure-flow relation for steady flow through the channel. This equilibrium relation is well characterized by the Forchheimer equation used to describe steady flow through porous media, where inertial effects become substantial when Re is of order 10 (34). The Forchheimer equation takes the form ΔP = αμLV + βρLV2, where V is the average velocity, α is the reciprocal permeability, and β is the non-Darcy flow coefficient. The latter two parameters are solely dependent on the system geometry and not on the working fluid. For our two-dimensional channel with obstacles, we take V = Q4/w so that the pressure-flow relation for the channel becomes P4=αμL4Q4/w+βρL4Q42/w2, where α and β are fit from simulations (Supplementary Materials and fig. S4). We account for this nonlinearity in a dynamical equation for Q4 by using a flow rate–dependent function in place of the constant resistance in Eq. 1. Specifically, we take R4(Q4) = αμL4/w + βρL4Q4/w2 to recover the Forchheimer equation in steady flow. A consequence of the nonlinearity of this channel is that it gives rise to a nonmonotonic relation between the pressure difference across the chicane channel and Pin. As Pin is increased from zero, Q3 initially increases, before decreasing, as indicated by the low-flow solution branch in Fig. 4B.

We now construct the dynamical model for the full network in Fig. 1 as follows: (i) we use flow relations of the form in Eq. 1 with constant resistances for the three channel segments without obstacles and with a flow rate–dependent resistance function (discussed above) for the obstacle-laden channel; (ii) we use Eqs. 2 and 3 to describe the flow rate and vortex dynamics in the chicane channel with ΔP34 substituted by (ζP3P4), where ζ is a free parameter that may deviate from 1 to account for an effective pressure difference across the chicane channel due to the finite size of the channel junctions; and (iii) we account for the most dominant minor pressure losses due to diverging flows at the channel junctions (35). For the latter, we include terms of the form kQ3Q5/Q1 in the flow equations for Q3 and Q5, where k is a positive constant. This leads to six ordinary differential equations (five for flow rates and one for the vortex size), which can be reduced to four equations by making use of the equations that account for flow rate conservation at the channel junctions: Q1 = Q3 + Q5 and Q2 = Q4Q3 (see Supplementary Materials for details of the model).

The model predictions of the total flow rate, chicane channel flow rate, and vortex size for the network in Fig. 1 under a common driving pressure at the inlets are presented in fig. S5. The model captures well the complex solution structure observed in Figs. 2A and 4 (B and C), shows strong quantitative agreement with simulations, and provides several interpretations for the observed flow behavior. First, spontaneous oscillations are found to arise through the transition from a fixed-point solution to a stable limit cycle via a supercritical Hopf bifurcation. The amplitude of the limit cycle grows with the driving pressure and eventually collides with the unstable solution surface of Q3 and r, as shown in fig. S5, thereby destabilizing the oscillating solution through a homoclinic bifurcation. Second, the nonlinearity arising from the Forchheimer effect gives rise to the two distinct bistable regions (and thus two negative-conductance transitions), as a result of the nonmonotonic relation between Pin and the pressure loss along the chicane channel. Third, the difference in the total flow rate between the solution branches is primarily determined by the minor losses. Without these terms, the model may still predict bistability, but the difference in total flow rate for solutions belonging to different branches would be negligible.

Our model can also be used to integrate the nonlinear behaviors described above into larger microfluidic systems. As an example, we consider an extended network with three outlets (Supplementary Materials and fig. S6), with two separate inlet flows. By driving the flows through this network using a common pressure, three unique oscillatory flow compositions can be realized at the outlets. Our model predictions show that the flow composition at the individual outlets is different, but the flow rate at the outlets oscillates in phase (Supplementary Materials and fig. S6B). Thus, the property that flow rates through all channel segments oscillate with the same period can be extended to larger networks and used to produce synchronized, time-dependent output flow patterns.


Motivated by the challenge of developing built-in controls in microfluidics, we identified mechanisms that can facilitate integration without dependence on movable parts or external actuation (other than through the working flow). This includes our demonstration of self-sustained oscillations, which can be used for timing and synchronization of flows through different channels; multistability and associated transitions, which can be used for signal amplification and switching; and hysteresis, which could serve as a possible mechanism for memory. In particular, we demonstrated the emergence of spontaneous periodic variations in the relative uptake rates from different inlets, which can be explored to generate time-dependent mixtures and output flow patterns. While these dynamical behaviors may resemble those found in microelectronics, they rely on effects that do not have direct analogs in electrical networks, namely, fluid inertia and the resulting nonlinearity arising from interactions between components.

Our results demonstrate that fluid inertia effects can be amplified and induce behaviors in fixed-structure microfluidic systems that have not been previously generated without external actuation. Indeed, the negative-conductance transitions, spontaneous oscillations, hysteresis, and multistability simulated and modeled in our system all emerge as a consequence of coupling between the geometric structure of the network and fluid inertia effects. Flows around obstacles and through the porous-like channel are determinant for generating these dynamics. Porous media microfluidics have become important for the study of flows through natural systems and laboratory-controlled experiments (36). In this work, we placed new emphasis on the viability of porous-like structures to serve as nonlinear fluid resistors and harness fluid inertia effects for nonlocal flow control throughout the network, which, crucially, can be realized as a built-in mechanism. Given that our system can be constructed from rigid materials, it is able to withstand a wide range of driving pressures (e.g., 1 to 106 Pa), which facilitates implementation across the length scales relevant to microfluidics.

The flow dynamics that arise in our system can be tailored for various applications. Microfluidic systems capable of carrying out sequential operations generally require a timing mechanism that is generated either from an external device or through the use of flexible valves (25). The oscillations that arise in our system could serve as an on-chip frequency reference and enable process synchronization or waveform synthesis. Moreover, the vortex dynamics that give rise to the oscillations may be used to enhance state-of-the-art methods for particle sorting and manipulation that function through interactions between particles and micro-vortices (37). In particular, vortex dynamics can be used to produce complex (and even chaotic) particle trajectories in laminar flows (38). Lastly, microfluidic networks are now widely used in the study of colloids (39) and active matter (40). Our system offers a rich environment to further investigate these materials given that they exhibit unusual collective behavior when placed in different flow fields (41) and when driven through porous media (42). Moving forward, we anticipate that the coupling of fluid inertia effects and network geometry can be further explored across microfluidic applications to yield new built-in flow control functionality.


Navier-Stokes simulations

Our simulations of the Navier-Stokes equations for incompressible fluid were performed using OpenFOAM-version 4.1. Meshes of the system geometry were generated using Gmsh version 2.9.3, with average cell area ranging from 10 to 70 μm2. The pisoFoam and simpleFoam solvers in OpenFOAM were used for time-dependent and steady-state simulations, respectively. For simulations where a Dirichlet static pressure boundary condition was used at an inlet/outlet, a Neumann boundary condition was used to set the gradient of the velocity field to zero in the direction normal to the inlet/outlet. This combination of boundary conditions results in a fully developed velocity profile at the inlet/outlet and corresponds to the physical situation in which the channels extend upstream and downstream of the computational domain. Similarly, when instead the flow rate was controlled at an inlet, a parabolic velocity profile was specified and a zero-gradient boundary condition was used for the pressure. Values for P3 and P4 in Fig. 5 were measured by averaging the pressure sampled across the channel width at a distance 3w/5 downstream of the chicane channel junctions.

Network dimensions

For the network presented in Fig. 1, the individual channel segment lengths, as labeled in Fig. 1D, are L1 = 0.1 cm, L2 = 0.6 cm, L3 = 0.1 cm, L4 = 1.0 cm, and L5 = 0.5 cm. The cylindrical obstacles in the obstacle-laden channel (Fig. 1C) are separated by a distance of approximately 6w/5, where w = 500 μm is the channel width (common to all channels in the network). The blade-like barriers in the chicane channel (Fig. 1B) are each placed a distance w/2 from the midpoint of the axis along the channel.

Reynolds numbers

The characteristic length scale used in defining the Reynolds numbers of the flows is the hydraulic diameter of the channels, defined as 4A/P, where A is the area and P is the perimeter of the channel cross section (common to all channel segments). In two dimensions, the hydraulic diameter is 2w, and the characteristic velocity used is Q/w. Therefore, we define the Reynolds number for individual channel segments to be 2ρQ/μ, where Q is the associated flow rate through the channel.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: Funding: This work was supported by ARO grant no. W911NF-15-1-0272 and a Northwestern University Presidential Fellowship. Author contributions: D.J.C., J.-R.A., and A.E.M. designed the research. D.J.C. performed the research. D.J.C. and A.E.M. wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from A.E.M.

Stay Connected to Science Advances

Navigate This Article