Report

How boundaries shape chemical delivery in microfluidics

See allHide authors and affiliations

Science  09 Dec 2016:
Vol. 354, Issue 6317, pp. 1252-1256
DOI: 10.1126/science.aag0532

Aspects of the design

In microfluidics systems, the small size of the channels ensures that the flow profiles are laminar, so solute mixing is largely governed by diffusion in the absence of active mixing. Aminian et al. revisited the classic phenomenon of Taylor diffusion to investigate the effects that the aspect ratio of the conduit has on the longtime axial distribution of solutes. They show both numerically and experimentally that the aspect ratio controls the skewness of this distribution, and thus pipe design alone is enough to control mixing profiles.

Science, this issue p. 1252

Abstract

Many microfluidic systems—including chemical reaction, sample analysis, separation, chemotaxis, and drug development and injection—require control and precision of solute transport. Although concentration levels are easily specified at injection, pressure-driven transport through channels is known to spread the initial distribution, resulting in reduced concentrations downstream. Here we document an unexpected phenomenon: The channel’s cross-sectional aspect ratio alone can control the shape of the concentration profile along the channel length. Thin channels (aspect ratio << 1) deliver solutes arriving with sharp fronts and tapering tails, whereas thick channels (aspect ratio ~ 1) produce the opposite effect. This occurs for rectangular and elliptical pipes, independent of initial distributions. Thus, it is possible to deliver solute with prescribed distributions, ranging from gradual buildup to sudden delivery, based only on the channel dimensions.

By exploiting the flexibility afforded by working at small scales, microfluidics and lab-on-a-chip devices have attracted interest because of their potential for simplifying and reducing costs of fluid-based laboratory processes. A key component of these devices is the pressure-driven transport of solute through channels. Recent technological advances, such as in microfluidic flow injection analysis (1, 2) and chromatographic separation (3, 4), have increased the demand for improved precision in the controlled delivery of solutes in such flows. This becomes particularly important when limited sample sizes are available, as in some applications envisioned for microfluidic laboratory analyses. We have found that the channel’s aspect ratio offers a method to specify the longitudinal solute distribution, including its peak locations and asymmetries: On long time scales, thin channels deliver species with a sharp front followed by a long tapering tail (“front-loaded”), whereas thick cross sections produce the reverse result (“back-loaded”) (see Fig. 1, A and B, for schematic illustrations).

Fig. 1 Setup, schematics, and distribution evolution.

(A) Flow geometry in pipes and ducts. (B) Back-loaded (left) versus front-loaded (right) distributions. (C) Snapshots of the Monte Carlo tracer distribution (T) projected through the longest cross-sectional direction (z) and (inset) the complete cross-sectionally averaged C(x,τ) for the circular pipe (left group) and rectangle λ = 0.2 (right group) with Pe = 104 at four times. Dashed lines denote the mean velocity; light blue lines mark the median of the distribution. The horizontal axis is scaled to the leftmost and rightmost particle at each τ, and peak intensity is set by the instantaneous scalar maximum, with colors indicating concentration from zero (black) < purple < red < maximum (yellow). Skewness values (color) and time (white) are indicated at the bottom and top of the insets. We observe more mass to the left (or right) of the mean in the circle (or rectangle) at later times.

Such phenomena go beyond the enhancement over molecular diffusion that was first explained and characterized by Taylor (5) (and has since been referred to as “Taylor dispersion”), with the prediction that flowing solute spreads with an effective, boosted diffusivity inversely proportional to the molecular diffusivity. This magnification is observed beyond long, diffusive time scales once the solute has diffused across the channel.

Recent research focusing on microfluidic applications has examined how the enhanced diffusivity also depends on the cross-sectional shape of the channel (69). An important yet largely unexplored avenue of investigation lies in how the cross section itself induces an asymmetry in the tracer distribution. For the special case of a circular pipe, it was observed theoretically that the distribution could arrive asymmetrically (10, 11), but the precise role played by geometry in establishing the symmetry features, particularly regarding the role that the aspect ratio plays in setting front- versus back-loading, appears to be unexplored in the literature.

On the very shortest time scales, the symmetry properties of solute distributions are nonuniversal with respect to the cross-sectional geometry (12): Elliptical domains preserve initial symmetries, whereas rectangular ducts instantaneously create asymmetries. In contrast, we show here that on longer time scales, the longitudinal asymmetries are universal with respect to a broad class of geometries and are essentially set only by the aspect ratio. We are primarily interested in channels with elliptical and rectangular cross sections (“pipes” and “ducts”), whose minor and major axes (i.e., short and long sides) are 2a and 2b, respectively, with the aspect ratio λ = a/b (0 < λ ≤ 1). The Péclet number, a nondimensional parameter measuring the relative importance of transport via flow versus diffusion, is defined on the basis of the shortest length scale, Pe = Ua/κ, and the diffusive time scale, td = a2/κ, where U is a characteristic flow speed and κ is the molecular diffusivity of the solute. Péclet numbers for microfluidic applications are typically between 10 and 105, with diffusivities ranging from 10–7 to 10–5 cm2/s (13). Thus, for these applications, typical length scales and flow speeds place relevant observations well beyond the nonuniversal behavior seen on the shortest time scales and onto intermediate-to-long time scales with respect to td.

After appropriate rescaling of variables (with physical time t normalized by td and denoted by τ, with longitudinal coordinate x, short transverse coordinate y, and long transverse coordinate z all normalized by a), the solute distribution T is governed by an advection-diffusion equation of the formEmbedded Image(1)with insulating (no-flux) boundary conditions at the wall and an initial condition depending only on the unbounded longitudinal (along channel length) variable x, which is taken to be a Dirac delta function, unless otherwise noted. The fluid flow u(y, z) is a laminar, steady-state solution to the Navier-Stokes equations with no-slip (no flow at wall) boundary conditions, driven by a negative pressure gradient. For the case of straight channels of fixed cross section, the fluid flow is expressible exactly in terms of either a polynomial or a Fourier series (see supplementary text for complete details). The partial differential equation in Eq. 1 couples the physics of fluid advection (represented by the second term) and molecular diffusion (the last term), whose solution at each time, τ, gives the spatial distribution of the concentration within the channel. Figure 1A depicts the setup and labels the various dimensions and coordinates.

A convenient tool for measuring asymmetry is statistical skewness, defined to be the centered, normalized third moment Sk(y, z, τ) = ∫(x – μ)3T dx/[∫(x – μ)2T dx]3/2, where μ = ∫xT dx is the mean of the distribution. This is the lowest-order integral statistic that measures asymmetry of a given distribution and whose sign typically indicates its front-loading (negative) or back-loading (positive).

Constructing this statistic requires the first three longitudinal moments of the distribution. The “partial” (pointwise) moments of the distribution Tn(y, z, τ) = ∫xnT(x, y, z, τ) dx are longitudinal moments at a given position within the cross section. These moments obey a similar set of equations to Eq. 1, called the Aris moment equations (14). Our main results depend on calculating long time behavior of the first three moments, which we can use to construct the skewness.

Figure 1C illustrates the particle distribution T computed by a Monte Carlo simulation for a circular and thin rectangular cross section for time increasing from short, advective-dominated times to longer, diffusive time scales. Shown in each panel’s inset is the cross-sectionally averaged distribution, C(x, τ) = |Ω|–1Ω T(x, y, z, τ) dy dz as a function of x, where |Ω| is the area of the cross-sectional domain. Initially, the circular pipe has approximately zero cross-sectionally averaged skewness (henceforth referred to as average skewness), becoming positively skewed at intermediate times and finally symmetrizing on long times. Alternatively, the particle distributions in the thin duct develop negative average skewness on short times and remain negatively skewed for all times. The strong relationship between the mean, median, and skewness further documents the notion that negative skewness implies front “loadedness,” as indicated by the relative position of the median and the mean. The horizontal axis is scaled by the maximum particle separation in the x direction (which roughly grows as the maximum flow speed times τ), capturing the short time spreading, whereas the distribution’s standard deviation grows slower Embedded Image at long time scales (15), which appears as a narrowing of the distribution at later times. This effect proves useful to highlight the complex interplay between short and long time regimes. Lastly, we observe that on longer time scales the circular pipe has a back-loaded distribution, whereas the thin duct is front-loaded.

We next compared these observations to laboratory experiments, performed as follows: First, we injected a small amount of fluorescein dye into a 300-mm-long square duct (1-mm–by–1-mm internal cross section) filled with distilled water, which was left to diffuse for several minutes to homogenize the initial condition. A syringe pump then pushed water through the duct, advecting the dye downstream. Side-view images were acquired with a digital still camera and postprocessed to construct the cross-sectionally averaged concentration as a function of downstream position. This procedure was repeated, matching the average flow velocity of the first experiment, on a thin duct (1 mm by 10 mm) (see supplementary materials for full details).

A direct comparison between Monte Carlo simulations and our laboratory experimental results performed in square and thin rectangular ducts is presented in Fig. 2, A and B, for intermediate time scales. Figure 2A documents the pointwise agreement between the experiment and the simulation viewed orthogonally to the long side. The near-uniform transport away from the top and bottom walls is also evident in both the experiment and the simulation; this is a consequence of flow in thin rectangles. This feature is absent for flow in elliptical cross sections and may be particularly useful in microfluidic channel design to reduce overall spreading. Figure 2B presents a comparison between the Monte Carlo simulation with the lab experiment for the cross-sectionally averaged tracer distribution at several different output times, illustrating that the thin domain is heavily front-loaded, whereas the square domain is strongly back-loaded (with variations among runs of typically less than 5%). The experimental observations quantitatively match the theory. The front- versus back-loading properties emerge as predicted by the theory in all variations of the experiment, including changes in initial data, thin-duct orientation with respect to gravity, flow rates, and concentrations.

Fig. 2 Experimental, computational, and theoretical comparisons.

(A) Thin-duct comparison: experiment (left) and simulation (right), projected down the short y axis (see Fig. 1A schematic), at τ = 0.2. (B) Experimental and computational comparison of the cross-sectionally averaged tracer distribution along the duct (top: square duct, cross section 1 mm by 1 mm; bottom: thin rectangular duct, cross section 1 mm by 10 mm) at τ = 0.2 and τ = 0.4. Average flow velocities are matched (square: Pe = 1670; thin: Pe = 752). Initial conditions were chosen to match the experiment (cross-sectionally uniform Gaussian in x, standard deviation: square, 3 mm; thin, 7.5 mm). (C) Skewness evolution: ellipses (left) and rectangles (right) in linear and log (inset) values, for varying aspect ratios. Equation 2 theoretical lines (all dashed) are shown for elliptical pipes and infinite (λ = 0) parallel plates. Special elliptical and rectangular cross sections of aspect ratio λ = λ* separate positive and negative relaxation with faster decay (dashed yellow lines).

The physical mechanism responsible for the net loading is a competition between the wall and interior. The shear flow shapes the concentration into a paraboloid, which creates an asymmetry: Away from the centerline, there is solute behind the leading front, providing a diffusive source that feeds the interior distribution and produces a left heavy tail (locally). This effect reverses at the wall. The relative influence of the wall versus the interior sets the net loading, and mathematical analysis is required to determine the precise balance.

Toward this end, we analytically assessed this dependence on the aspect ratio through the calculation of the first three moments. Each moment is an evolving function of the location within the cross section, satisfying a driven heat equation with a driver that depends on the fluid flow and lower moments. Each moment’s long time behavior is reduced, for its cross-sectional average, to a simple ordinary differential equation and a domain-dependent partial differential equation for its fluctuation (see supplementary text for complete details).

These calculations yield a simple formula for the long time skewness as a function of the aspect ratio and specific domainEmbedded Image(2)where the functions G1 and G2 are defined in the supplementary text and depend only on the aspect ratio and domain shape and not on location within the cross section. This surprising result predicts that the long time nonzero skewness becomes independent of y and z. The coefficient of the decay rate τ–1/2 can be calculated explicitly in closed form for elliptical cross sections and agrees with previous exact results for the circular pipe (10, 11).

The results of Monte Carlo simulations are shown in Fig. 2C, with varying aspect ratios in elliptical and duct domains. The majority of plots show this τ–1/2 decay rate at long times with a few exceptions. For the very thin domains (λ = 0.01), the long time regime is not yet reached, as the final diffusive time scale is Embedded Image or τ = 104. Notably, there are critical aspect ratios for the ellipse and duct domains that exhibit faster decay and separate positive (back-loaded, wall-dominated) from negative (front-loaded, interior-dominated) relaxation in time.

Our theory can directly predict these “golden,” cross-over aspect ratios: Because the denominator of Eq. 2 can be shown to be strictly positive, the long time sign of the skewness is determined by the numerator G2.

To handle the general elliptical case, working in elliptical coordinates yields exact solutions for G1 and G2, in terms of a finite Fourier sum (see supplementary text for full calculations)Embedded Image(3a)Embedded Image(3b)[The expression for G1 is proportional to the coefficient for classical diffusive enhancement and agrees with previously calculated exact results for the ellipse (14, 16)].

The second term G2 provides an exact prediction for the cross-over (“golden”) aspect ratio, separating front-loaded from back-loaded behavior at long times. The expression for G2 has a simple, irrational root at aspect ratio Embedded Image. At this value, the skewness vanishes at a faster rate ~τ–3/2 instead of ~τ–1/2, which we observe in the numerics (Fig. 2C) with a reference –3/2 slope line. Further, the ratio 3G2/(2G1)3/2 provides the correct large Péclet prefactor to the long time decay rate as a function of aspect ratio for the skewness (defined in Eq. 2) shown by the predicted reference lines in the left inset of Fig. 2C. The agreement between the simulation and these theoretical curves can be seen from Fig. 2C to hold for time scales well before the Taylor regime is reached. In general, the long-time leading order asymptotics will be ~Aτ–1/2 + Bτ–3/2, where A is given in Eq. 2 and B is a known function of the aspect ratio and Péclet number. For aspect ratios near but different from the golden ratio, A will be very small and the scaling τ–3/2 will be observed for a large period of time.

As with the family of ellipses, there is a distinguished “golden” rectangle. For this, we employ a finite element scheme to solve the partial differential equations governing the moments at long times and calculate G2 as a function of the aspect ratio. The distinguished rectangle aspect ratio is Embedded Image. Again, at this value, the skewness decays as ~τ–3/2, which agrees with simulations shown in Fig. 2C. Note the success of the theory in quantitatively predicting the skewness well before the final diffusion time scale. Aspect ratios below Embedded Image have negative long time skewness (front-loaded), and those above Embedded Image are positive (back-loaded).

Theoretical results such as these are useful in designing optimal microfluidic devices. We emphasize that our analysis gives an exact result for the elliptical domains and an accurate numerical prediction in the rectangular duct. The fact that these critical aspect ratios are so close for quite different domains further emphasizes the universality of the theory, which is particularly important in experimental applications where geometric variations are inevitable.

The role of geometry in front- and back-loading can be visualized computationally by observing the evolution of the cross-sectionally averaged distribution C(x,τ). Figure 3 shows this evolution for varying aspect ratios (thin to thick, left to right) and geometries (top: ellipse; bottom: duct), using stacked snapshots of the distribution with time increasing from top to bottom. Again, thin domains produce front-loaded distributions, whereas thick geometries reverse this result. Also note the difference between elliptical and rectangular domains at short and intermediate times, in which the ellipse experiences symmetric spreading in contrast with the rectangular case.

Fig. 3 Evolution of the distribution C(x,τ).

Ellipses (top) and rectangles (bottom) for aspect ratios of λ = 0.2, λ*, and 1 illustrate nontrivial asymmetry depending on geometry over many time scales (time flows from top to bottom). Distributions are scaled by their peak value and standard deviation for visualization. The mean velocity is marked by vertical lines. The front-loading evolution can be seen in the left column, whereas the back-loading persistence is apparent in the right column. The middle column (golden ratios) illustrates more rapid symmetrization.

Although difficult to probe in a lab experiment, we can computationally explore pointwise distributional variations within the cross section itself, as opposed to its average. Figure 4 presents the evolution of pointwise skewness as a function of aspect ratio and geometry. In agreement with our general theory, at long times, the skewness evolves to a uniform state throughout the cross section. On earlier times, there is interesting dynamics connecting the short-time positive skewness (red) at the walls with negative skewness (blue) in the interior and subsequent evolution to the long time state. This demonstrates the complex competition between the wall and the interior, and the interior ultimately prevails in the thin domains shown here.

Fig. 4 Pointwise skewness over cross sections.

Parameters are Pe = 104 and aspect ratios (left to right) of λ = 0.2, λ*, and 1.0 at times τ = 0.0014, 0.008, 0.046, 0.44, and 2.5 (top to bottom). Short time skewness: positive at the walls (red) and negative in the center (blue) in all cases due to diffusive pumping. Pointwise skewness for all cross sections is approximately uniform at τ = 2.5 (compare with Eq. 2), with golden ratios λ* separating positive from negative values.

The front- versus back-loading effects we have documented are expected to be an important factor and will ultimately influence the design of microfluidic devices by affording control of mass distribution by simple geometric means. This result may be particularly useful in flow injection analysis, where a sample bolus is injected into a carrier fluid that reacts and generates by-product gradients that are analyzed downstream (1, 2). Our new geometric control of concentration profiles allows for the tuning of gradients to optimize output response. This effect could further be used to control and generate specific chemical gradients within and at the exit of microfluidic channels (17), which may be useful for chemotaxis (1820) and drug development (21) assays. Moreover, the additional control could improve applications such as open-channel or hydrodynamic chromatographic separation processes (3, 4), where longitudinal asymmetry presents considerable measurement challenges (22, 23). For instance, the aspect ratio phenomenon we have discovered could be exploited to minimize or actively reverse the asymmetries encountered during the separation process, sometimes referred to as “peak fronting” or “tailing” (23, 24). Among the various yet-unanswered questions are whether more complex geometries can provide additional control and, in particular, whether more extreme asymmetries can be obtained.

Supplementary Materials

www.sciencemag.org/content/354/6317/1252/suppl/DC1

Materials and Methods

Supplementary Text

Fig. S1

References (25, 26)

References and Notes

Acknowledgments: We acknowledge funding from the Office of Naval Research (grant DURIP N00014-12-1-0749) and the NSF (grants RTG DMS-0943851, CMG ARC-1025523, DMS-1009750, and DMS-1517879). The experimental methodology and the details of the Monte Carlo simulations are documented in the materials and methods and the supplementary text.
View Abstract

Navigate This Article