## Abstract

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.

## INTRODUCTION

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 (*4*–*6*). 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 (*9*–*11*). 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 (*15*–*18*), 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 *P*^{out} to be held at a fixed common value, taken to be zero. At the inlets, we control either the pressures (*Q*_{1} and *Q*_{2}).

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/m^{3} 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,

In Fig. 2A, we show simulation results of the total flow rate *Q _{T}* =

*Q*

_{1}+

*Q*

_{2}through the network in Fig. 1 over a range of driving pressures,

*P*

^{in}, from which we observe two notable properties. First, for

*P*

^{in}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

*P*

^{in}being fixed. In particular, we find that at a critical value of

*P*

^{in}, solutions along the high-flow branch (red symbols in Fig. 2A) become small-amplitude limit cycles. The corresponding amplitudes and periods grow with

*P*

^{in}(the frequency of the oscillations decreases from 20 to 4 Hz; see the Supplementary Materials and fig. S2). At a higher critical

*P*

^{in}, 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.

Another outstanding property that arises from the bistability in our system is the possibility of negative-conductance transitions and other sudden transitions in *Q _{T}* that result from small changes in

*P*

^{in}. We characterize these transitions, which occur at the boundaries of the bistable regions (Fig. 2A), by defining (local) fluidic conductance and resistance as

*C*= δ

*Q*/δ

_{T}*P*

^{in}and its reciprocal, respectively. Here, δ indicates a finite change and

*P*

^{in}is the controlled variable. Therefore, negative-conductance and negative-resistance transitions occur when an increase (decrease) in

*P*

^{in}leads to a decrease (increase) in

*Q*. Notably, as shown in Fig. 3A, our system exhibits transition points at which

_{T}*C*(δ

*P*

^{in}) diverges in the limit of small δ

*P*

^{in}: two points at which

*C*(δ

*P*

^{in}→ 0) = +∞, corresponding to positive-conductance transitions, and two points at which

*C*(δ

*P*

^{in}→ 0) = −∞, corresponding to negative-conductance transitions. Figure 3B shows that related transitions emerge when the flow rate

*Q*, rather than the pressure

_{T}*P*

^{in}, is taken as the control variable. In this case, a change in

*Q*can lead to transitions in which

_{T}*P*

^{in}changes by a finite amount. In particular, the latter includes signal amplification transitions, which are remarkable transitions in which an infinitesimal increase (decrease) in

*Q*leads to a finite decrease (increase) in

_{T}*P*

^{in}. 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.

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 *Q*_{3}. 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 *Q*_{3} 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 *Q*_{3} and *r*.

We determine the relationship between *r* and *Q*_{3} by performing simulations in which the flow rates at both inlets (*Q*_{1} and *Q*_{2}) are controlled. From these simulations, we compute *r*, *Q*_{3}, and the pressure loss along the chicane channel Δ*P*_{34}, where the latter corresponds approximately to *P*_{3} − *P*_{4} (Fig. 1). In Fig. 5, we show relations between these quantities for sets of simulations in which *Q*_{1} is fixed while *Q*_{2} is varied. We observe nonlinear relations between *r* and *Q*_{3} (Fig. 5A), between *Q*_{3} 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 *Q*_{2} 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.

### 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 as*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*/*w*^{3} 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, *R*_{3}, and the vortex size, *r*, for 60 < *r* < 400 μm (Fig. 5C). To construct an approximate dynamical equation for *Q*_{3}, we use the form of Eq. 1 with the constant resistance replaced by a function of *r*. Specifically, we take *R*_{3}(*r*) = 12μ(*L _{b}* + γ(

*r*−

*r*)

_{b}^{2})/

*w*

^{3}, where

*L*serves as a base component of the resistance, γ is a constant coefficient of the variable component that depends on the vortex size, and

_{b}*r*is the vortex size that minimizes the resistance (from Fig. 5C,

_{b}*r*≈ 150 μm). With this added dependence on

_{b}*r*, we must also account for the dynamics of the vortex size. The steady-state relation between

*Q*

_{3}and

*r*found through flow-controlled simulations (Fig. 5A) can be well fit by a cubic equation of the form

*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 form

*P*

_{34}, 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

*Q*

_{1}can be accounted for by allowing η and

*L*to be functions of

_{b}*Q*

_{1}(Supplementary Materials).

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* + βρ*LV*^{2}, 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* = *Q*_{4}/*w* so that the pressure-flow relation for the channel becomes *Q*_{4} by using a flow rate–dependent function in place of the constant resistance in Eq. 1. Specifically, we take *R*_{4}(*Q*_{4}) = αμ*L*_{4}/*w* + βρ*L*_{4}*Q*_{4}/*w*^{2} 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 *P*^{in}. As *P*^{in} is increased from zero, *Q*_{3} 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 Δ*P*_{34} substituted by (ζ*P*_{3} − *P*_{4}), 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 *kQ*_{3}*Q*_{5}/*Q*_{1} in the flow equations for *Q*_{3} and *Q*_{5}, 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: *Q*_{1} = *Q*_{3} + *Q*_{5} and *Q*_{2} = *Q*_{4} − *Q*_{3} (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 *Q*_{3} 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 *P*^{in} 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.

## DISCUSSION

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 10^{6} 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.

## MATERIALS AND METHODS

### 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 μm^{2}. 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 *P*_{3} and *P*_{4} in Fig. 5 were measured by averaging the pressure sampled across the channel width at a distance 3*w*/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 *L*_{1} = 0.1 cm, *L*_{2} = 0.6 cm, *L*_{3} = 0.1 cm, *L*_{4} = 1.0 cm, and *L*_{5} = 0.5 cm. The cylindrical obstacles in the obstacle-laden channel (Fig. 1C) are separated by a distance of approximately 6*w*/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 4*A*/*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 2*w*, 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 MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/20/eaay6761/DC1

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.

## REFERENCES AND NOTES

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

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).