Chapter 4: Hele-Shaw Flow Theory in the Context of Open Microfluidics: From Dipoles to Quadrupoles – Open-Space Microfluidics

Chapter 4
Hele-Shaw Flow Theory in the Context of Open Microfluidics: From Dipoles to Quadrupoles

Étienne Boulais1 and Thomas Gervais1,2

1École Polytechnique de Montréal, Department of engineering physics, Montreal, Canada

2Centre de recherche du Centre hospitalier de l'Université de Montréal (CRCHUM), Montreal, Canada

4.1 Introduction

Henry Selby Hele-Shaw1 (1854–1941) was a British mechanical engineer and inventor who made major contributions to the fields of fluid dynamics, aeronautics, and automobile engineering. Some of his inventions include the Hele-Shaw clutch (1905), the first friction clutch for automobiles, and the variable pitch propeller (1924) for aircrafts. But perhaps his most important contribution to science and technology is an abstract device called a “Hele-Shaw cell,” initially a method that he invented to help his students visualize streamlines of a flow by confining fluids between two closely spaced plates (Figure 4.1). In a paper published in 1898, Hele-Shaw further demonstrated that, if the plates in his cell were close enough, flow was laminar at all velocities [1].

Figure 4.1 Dye injection and streamlines around a (thick) flat plate in a macroscopic Hele Shaw configuration with a gap of 1 mm. Flow from left to right, dye injection and streamlines around a (thick) flat plate, horizontal setup, gap ∼1 mm. The “plate” is modelled with a rubber pad. Photo courtesy of Hubert Chanson (University of Queensland).

Nowadays, Hele-Shaw flows are viewed more generally as a class of quasi-two-dimensional (2D) flows, arising from lubrication theory, and mathematically akin to Darcy flows in porous media and to inviscid flows [2]. Their most striking feature is that, contrary to the ordinary Stokes flow, they can be described by a single scalar function, or velocity potential, from which all flow properties can be calculated. In this sense, Hele-Shaw flow fields bear a strong mathematical analogy with electrostatic, magnetostatic, and gravitational fields. Thus, the same mathematical tools and principles that apply to these classical fields of physics can also be applied to study Hele-Shaw flows. There is a rich literature on advanced mathematical methods to solve for these types of quasi-irrotational flow fields using Green's functions and complex variable techniques such as conformal mapping [3] or multipolar expansion [4].

In this chapter, we will first describe the approximations under which the classical Navier–Stokes equation reduces to the Hele-Shaw problem (Section 4.2). We will also describe and discuss the validity of the point source approximation in Hele-Shaw cells and how they relate to experimental results. The limit and validity of these approximations will provide design and operation criteria to ensure that the probes will behave as predicted by the Hele-Shaw theory. Then, we will formulate the microfluidic dipole and quadrupole problems using Hele-Shaw flows (Section 4.3). The validity of the analytical results will be compared with full three-dimensional (3D) numerical simulations of the Navier–Stokes equation in the same domain and, to some extent, with already-published experimental data on probe operation. Next, we will add one layer of complexity and use the microfluidic probe (MFP) flow profiles to study diffusion at the edges of the probe's hydrodynamic flow confinement (HFC). Key approximations will be provided to extract scaling laws and accurate 1D models to quantify diffusive transport under MFPs (Section 4.4). Finally, we will conclude with a brief overview of how these models can be practically applied to optimize probe design, guide experiments, and automate probe operation (Section 4.5).

4.2 Fundamentals of Hele-Shaw Flows

In this section, we describe the underlying assumptions used to derive Hele-Shaw flows from the general Navier–Stokes equation. In doing so, we describe the limits of the Hele-Shaw flow approximation and illustrate the basic features of this kind of quasi-2D flow profile.

4.2.1 Derivation of Hele-Shaw Equation from the Navier–Stokes Equation

The Navier–Stokes equation is typically the entry point to virtually all branches of fluid mechanics. This nonlinear vectorial time-dependent partial differential equation (PDE) is well established, and its solutions, with the appropriate initial and boundary conditions, arguably describe every possible physical problem in laminar fluid dynamics. In the absence of external forces other than pressure, its form is well known [5] and given by


where ρ is the fluid density, μ the fluid's dynamic viscosity, p the pressure distribution, and the 3D vector field. Therefore, we shall demonstrate in this section a well-known result that Hele-Shaw flows are a subset of approximate solutions of the full Navier–Stokes problem obtained by exploiting simplified symmetries in quasi-2D geometries. First and foremost is that, in Hele-Shaw theory, flow occurs in between narrowly spaced infinite parallel plates. Using a Cartesian reference frame, it is therefore restricted in the vertical direction to the plate gap and will occur over an undetermined characteristic length L in the planar directions from an arbitrary origin O (, ).

To highlight the Hele-Shaw flow behavior, Navier–Stokes expression is better expressed in terms of dimensionless quantities using the following dimensionless groups: , , and , , , , . Expanding the N–S equation, we obtain three dimensionless scalar equations of the form

These equations greatly simplify when three key conditions are met:

  1. 1. The gap is much smaller than the characteristic length L over which the flow pattern of interest occurs (). This condition implies that the curves of the streamlines in the x and y direction are much smoother than those in the z direction. In this case, the velocity profile in z can always be assumed to be parabolic locally.
  2. 2. Reynolds number is much smaller than 1, such that we can approximate the flow as purely creeping (Stokes flow, ). Note that the characteristic length used in the Reynolds number is given by , also used for microfluidic channels [6]. This result implies that for flow patterns in between parallel plates occurring over a sufficiently long distance L, creeping flow can always be assumed.
  3. 3. There is no significant pressure difference occurring in the z direction (). In practice, this condition is always met in microfluidic systems as the pressure head from a fluid in a small gap G, , is negligible. Finally, the no flux boundary conditions in z at and , directly yields .

Using these three approximations, the terms written within braces in Eq. (4.2) are vanishingly small, and the Navier–Stokes equation simplifies to

In Eq. (4.3), the pressure is assumed not to vary in z and is thus a function of only x and y, p(x, y). The same can be said for the direction of the velocity vector, which can be written as . Since the pressure does not vary in z, integrating Eq. (4.3) twice over is easily done (using no slip boundary conditions for vx and vy at and ) and yields


This is the normalized Hele-Shaw flow expression. Denormalizing, we get the well-known expression of flow in a Hele-Shaw cell (where the origin is located at the bottom plate) [5]:

The two components of Eq. (4.5) can be easily combined to form a 2D velocity field of the form

The quasi-2D nature of Hele-Shaw flows becomes obvious under this form as

  • The z component of the velocity and of the pressure gradient is always zero.
  • The velocity vector can be expressed as a product of a function f(z) (the parabolic flow profile in z) and a function g(x, y) (the pressure gradient function in x, y).

This results in a vorticity vector with a zero component in the z direction:

Thus Hele-Shaw flows are often viewed as quasi-2D irrotational flows. As seen from Eq. (4.7), they are purely irrotational at the midplane . Another way often used to represent Hele-Shaw flows as purely 2D flows, and that will be used in the following sections, is to use a height average of the flow. Averaging the velocity profile and vorticity from Eqs. (4.6) and(4.7) yields two key results:

Height-averaged Hele-Shaw flows can now be considered irrotational and purely 2D. Furthermore, if the fluid studied is incompressible (such as water), conservation of mass in any control volume requires that the divergence of the flow field be zero () in the absence of source terms. Applying this condition to Eq. (4.8) yields a 2D differential equation of the Laplace form


When considering the presence of a distribution of sources and drains, the Laplace equation generalizes to a Poisson-type equation of the form

where q(x, y) is a function describing the 2D spatial distribution (surface density) of sources and drains in the plane.

This has profound impact on the way they can be treated mathematically. Eq. (4.10) is a PDE of the Poisson form, mathematically analogous to the relationship between electrostatic potential from a charge distribution, electric potential and current density in a conductor, and gravitational potential and its related mass distribution. Due to conservation of mass (or charge), these fields are irrotational. The field can also all be readily computed using Gauss' theorem for distributions with obvious symmetries. Table 4.1 summarizes the mathematical analogy between these four situations.

Table 4.1 Mathematical analogy between 2D electrostatics, electrokinetics, Hele-Shaw, and Darcy flows

Electrostatics (charge qe) Electrokinetics (current I) Hele-Shaw flows (flow rate Q) Darcy flows (flow rate Q)
Potential (scalar) Electrostatic potential φ(x, y) Electrostatic potential φ(x, y) Velocity potential φ(x, y)
(or pressure)
Velocity potential φ(x, y)
(or pressure)
Field (vector) Electric field
Current density

σ: electrical conductivity
Velocity profile
Velocity profile
Poisson equation
Gauss's theorem (continuity equation)

In all four situations presented in Table 4.1, the field lines around the point source in planar 2D geometry are of the form


where is the vector field (electric field, current density, velocity field), A is a proportionality constant (qe/ϵ0, I, Q), G is the domain thickness, is the distance from the origin, and is the position of the point source with respect to the origin. The variables x, y, x′, and y′ play the same role in Cartesian coordinates. The potential function associated with this type of conservative field is, dropping an integration constant,

Incidentally, the potential function in Eq. (4.12), since it is the response to an impulse function (a point source term), is also the Green's function associated with this type of conservative fields in infinite planar 2D space. Written in Cartesian coordinates,

We shall review, in the next section, how this function can be convolved with any source and sink distribution to compute the Hele-Shaw velocity profile of any Hele-Shaw cell, with emphasis on flow under microfluidic dipoles and quadrupoles.

4.2.2 Hele-Shaw Point Sources, Round Monopoles, and Square Monopoles

Point flow sources are of course never achieved in practice. Flow has to be brought into the Hele-Shaw cell, either from above or from below, through finite openings that can have any shape and that perturb the flow locally. Yet, the point source approximation is extremely useful to study flow under MFPs as we shall see in this section.

Green's function in Eq. (4.13) can be convolved with any 2D source distribution term q(x′, y′) to yield the resulting velocity potential in the cell using the convolution theorem

For the case of flow outside circular monopoles of radius a, the method above can be applied, but the velocity profile is more readily obtained using Gauss's theorem instead:

The result obtained in Eq. (4.15) is identical to the flow around an ideal point source obtained by taking the gradient of the impulse response in Eq. (4.14). In other words, when outside a circular monopole, the flow does not depend on the aperture size and is identical to that of a point source located at the center of the circular monopole.

In many practical circumstances, fluidic apertures used to inject or aspirate fluids from a Hele-Shaw cell are better made square or rectangular rather than circular. That is the case when a planar flow profile near the aperture is desired [7] or when microfabrication processes favor openings with sharp corners [8]. The flow profile around such openings can once again be computed using Eq. (4.14). In the common case of a square aperture of side a, the source distribution term is


The convolution integral then yields


An analytic expression exists for the above equation and is plotted for visualization in Figure 4.2.

Figure 4.2 Comparison between flow around a square and a round aperture in a Hele-Shaw cell. (a) Plot of the isopotential lines near a square aperture (solid line) and a round aperture (dotted lines) in normalized coordinates (2x/a, 2y/a). (b) Normalized velocity profiles for square and round apertures along a line taken perpendicularly to the side of the square and crossing the origin at the center of the square. Relative error . The error for using a point source approximation on a square opening is down to 1.5% at a distance of only from the center of the aperture. (c) Full 3D FEM simulation (COMSOL 5) comparing isopotential lines from a square opening (solid) with a round opening (dashed) of same characteristic dimension a. The difference between the two opening types is imperceptible outside of the square aperture.

A key observation is that, when sufficiently far away from a square opening, or an opening of any shape of characteristic length a, the far field solution decays quickly to asymptotically match that of a round aperture, which is in turn identical to that of a point source aperture. In the particular case of a square aperture, the error for using a point source model is less than 1.5%, only half a width away from the aperture (Figure 4.2b). This makes the point source approximation extremely relevant in the study of flow under MFPs, and this approximation will be abundantly used in the next sections.

4.3 Applications to Microfluidic Dipoles and Quadrupoles

In the previous section, we have laid the theoretical background that allows the definition of a Hele-Shaw cell and that shows the impulse response of a point source located inside the cell. The flow under any MFPs composed of any number of injection and aspiration apertures can now be modeled using the associated Green's function and a source distribution term. In this section, we will first derive the flow profiles typical under the most common probe configurations (dipoles and square quadrupoles) and discuss how key flow parameters, such as stagnation points, HFC area, pressure, and shear stress, can be systematically investigated using the Hele-Shaw approximation. We will also discuss under which practical conditions the simple model developed here fails and to what extent.

4.3.1 Velocity Potentials for Dipoles and Quadrupoles

MFPs, as described earlier in Chapters 1 and 2, are composed of a flat mesa (or head) of dimensions in which two or more apertures of size a are bored and through which fluid is either injected or aspirated. The apertures are separated by a characteristic distance d. In a microfluidic dipole, d is the interaperture distance (center to center). In a square microfluidic quadrupole, d is the diagonal interaperture distance (center to center). The probe is brought within a distance G from a flat surface (where ), thus forming a narrow gap between the probe mesa and the surface where the fluid can circulate. The mesa is chosen to be much larger than the HFC such that its edges do not perturb the HFC zone and, in this analysis, will be considered infinite. (In practice, this is achieved when )

Figure 4.3 describes the general problem of representing dipoles and quadrupoles with point sources in 2D geometries (viewed from above). The total velocity potential and consequently the resulting velocity profile are computed by summing the contribution of every source term and the point of interrogation (x,y).

Figure 4.3 Modeling a dipole (A–C) and a square quadrupole (D–F) with distributed point sources. (A,B): The 3D probe geometries modeled in this chapter. (C, D) Equivalent 2D Hele-Shaw model used for analysis (replacing finite apertures by point sources). (E, F) Typical simulations obtained using full 3D simulations of the advection–diffusion problem in steady state using COMSOL. Arrows highlight the parabolic flow profile in the z direction. Solid lines are flow streamlines. Grayscale represents reagent concentration (, ). The position of stagnation points is identified. (Simulation values: , , , ).

For the dipole and quadrupole presented in Figure 4.3, the source distribution is given by


Using these values in the convolution integral, Eq. (4.14), the velocity potential describing the flow under both types of probes, after a few simplifications, is found to be

where the parameter is the ratio of aspiration to injection flow rates. From these potentials, all key quantities defining flow under dipoles and quadrupoles can be obtained.

4.3.2 Deriving Key Operation Characteristics for Dipoles and Quadrupoles Stagnation Points and the Hydrodynamic Flow Confinement Zone

Perhaps the most important characteristic of MFPs, as defined in the seminal paper by Juncker et al. [9], is the concept of HFC. In any probe configuration, when the net mass balance under the probe is zero (the sum of the flow rates injected through all injection apertures equals the flow rate aspirated back), the flow reaches to infinity (or the edges of the probe in practice) and the injected fluids are not confined. Yet, if the aspiration flow slightly exceeds the injected flow, as described by the condition , a net flow rate is aspirated from the medium below the probe. In the far field regime (), the probe can be viewed as a single point drain of strength . However, close to the injection aperture, the outward flow velocity is high (tending to infinity as one approaches a point source). Therefore, there must exist at least a point in space for each aspiration aperture where the velocity transitions from outward to inward, which is a point of zero velocity or stagnation point. By symmetry, this point has to be located on the probe's axis, as any point away from this axis will have a net y-velocity, which makes them relatively easy to compute.

For the case of the dipole, the position of the dipoles can be described, as in Table 4.2, by computing the velocity profile on the probe's axis (the x-axis) and finding the point where it is identically zero. By symmetry, the y-component of the velocity will be zero everywhere on the x-axis:


Table 4.2 Relationship between probe design and control parameters and the velocity potential function

Property Operation on phi Scaling law
Velocity potential φ
Pressure p
Stream functions (streamlines)
Average velocity profile vave
3D Hele-Shaw flow profile
Shear stress τ along walls in the flow direction
Position of stagnation points xSPd

From which we obtain

Looking at Eq. (4.22), we observed that the position of the single stagnation point in the dipole configuration is indeed at infinity when the net flow rate is zero under the probe . As the flow rate ratio α increases, the stagnation point is brought closer to the injection aperture until, when α is sufficiently large, it disappears within the injection aperture of finite diameter a (occurs when ).

Similarly, for the quadrupole, taking the gradient of Eq. (4.20) along the line of symmetry, the x-component yields

Eq. (4.23) has three roots

corresponding to the three stagnation points, all along the x-axis, in the square quadrupole configuration (at the center and outside each injection apertures). Once again, for a symmetric quadrupole (), the stagnation points are found at infinity (Figure 4.4).

Figure 4.4 Streamlines under a microfluidic dipole probe (a–d) and a microfluidic quadrupole (e–h) for various flow rate ratio α. Full 3D simulations using COMSOL. (, , ).

The results obtained in Eqs. (4.22) and (4.24) reveal that the position of the stagnation points under a given MFP is not directly controlled by either the injection or the aspiration flow rates, but rather by a single variable: the ratio of aspiration to injection flow rates. Therefore, there exists an infinity of injection and aspiration flow rates that will yield the same HFC area (provided that their ratio remains constant). In practice, this feature makes variables such as pressure, velocity, and shear stress, which are driven by the magnitude of the flow rates, entirely independent from the shape of the probe. In other words, the probe operator can specify both the flow velocity under the probe and the HFC independently. Any HFC forms can be achieved at low shear stress. Alternatively, flow can always be made fast enough for a given HFC that diffusive blurring becomes negligible around the probe, as will be seen in the next section.

4.3.3 Numerical Investigation of Model Accuracy

In order to develop the simple Hele-Shaw models for MFPs described in this section, a certain number of assumptions had to be made, which we recall here: (i) the gap between the probe and the surface has to be much smaller than d, the probe interaperture distance (Hele-Shaw approximation #1, Section 4.2.1); (ii) aperture size also has to be much smaller than the interaperture distance d (point source approximation, Section 4.2.2); and (iii) there is no vertical flow (z flow) under the probe as plates are separated by a constant gap G everywhere (Hele-Shaw approximation #3, Section 4.2.1). The “creeping flow” (low Re number) approximation (Hele-Shaw approximation #2, Section 4.2.1) is assumed always respected; thus far, no papers on MFPs have reported otherwise in the literature. However, in practice, these assumptions are not always respected. The most striking example being the first MFP paper published [9], where the probe design uses 25 µm apertures separated by an interaperture distance of 50 µm and a probe-to-surface gap of 10 µm (). In these circumstances, the effect of non-negligible aperture size should be felt under the probe. Similarly, the “no vertical flow” condition is not respected in the near vicinity of the aperture (or right under it) as streamlines comes from the probe surface and have to curve to enter the gap under the probe. Finally, the aperture holes in the probe create directly under the apertures a gap much larger than G and thus influence the streamlines and general hydrodynamic flow resistance from the inlet to the outlet. These effects of finite gap and aperture size have been singled out and analyzed previously [10]. Results compiled from 20 different full 3D numerical simulations and 3 experimental conditions, as developed previously [11], are also compiled in Figure 4.5. They show that the point source approximation slowly breaks down as . The main physical reason for this breakdown is the fact that, underneath flow apertures, flow is neither uniformly distributed under the flow surface, and hydraulic flow resistance is much smaller than inside the gap. Therefore, the length traveled in straight line from the inlet to the outlet inside the gap will be shorter by a distance of two aperture radii, or a, to effectively become (da) instead of d. This will have the effect of favoring flow along this straight path of lesser resistance, with the consequence of bringing the outer stagnation points closer to the inlets than expected. Square openings will also induce a greater perturbation that round openings as the ratio a/d increases. Finally, for the most outlying point in Figure 4.5 (for a square aperture with ), the stagnation point is effectively located under the aperture, where the Hele-Shaw theory breaks completely. Yet, in general, the position of stagnation points in Figure 4.5 is in close agreement with the Hele-Shaw model when a/d remains small (i.e., ).

Figure 4.5 Comparison between analytical models using Hele-Shaw theory and full 3D simulations of convective transport using the Navier–Stokes equation. In gray: normalized stagnation point positions as a function of flow rate ratio for microfluidic dipole probes with different aperture size ratios. Theory (solid line) compared with 20 different numerical results (markers). , . In black: normalized outer stagnation point positions as a function of flow rate ratio in microfluidic quadrupole probes. Theory (solid line) compared with experimental data (markers). , , .

4.4 Diffusion in Hele-Shaw Flows

While convection dictates the general shape of the flow profile and controls important parameters to MFP operations such as flow velocity, pressure, shear stress, and general HFC shape, they are insufficient to fully characterize MFP behavior. MFPs, as most microfluidics-in-the-open technologies, are designed to deliver reagents near surfaces, may they be fluorescent markers, antibodies, enzymes (such as trypsin), or stains (such as hematoxylin and eosin (H&E)). Typically, the reagents delivered are small molecules exhibiting significant Brownian motion. In such cases, reagent diffusion along the edges of the HFC will be non-negligible and will induce a “blurring effect” at its boundaries. This effect can be either highly nondesirable, when trying to precisely control reagent patterning onto surfaces, for example, or highly desirable to create precise, stable, and sharp concentration gradients. The amplitude of diffusion, or blurring effect, under an MFP is more challenging to model than the flow rates described in Section 4.3 as it depends both on convective and diffusive effects. Intuitively, if the flow under the MFP is made slow enough, the MFP should act like a source of constant concentration and reagents should slowly diffuse outward toward infinity. On the other hand, if the flow rate is made fast enough, gradients will be infinitely sharp at the probe's edges and diffusion will be considered negligible. This behavior, as we shall see, is typical of transport processes controlled entirely by the flow's Péclet number (the ratio of the diffusive to convective time scales) [12]. In this section, we will formulate the advection–diffusion transport equations under an MFP and highlight useful flow regimes to either exploit diffusion broadening for gradient generation or suppress it for precise surface patterning.

4.4.1 Advection–Diffusion Transport Equations

Solving advection–diffusion problems, both analytically and numerically, is a two-step problem, as it requires first that the streamlines, or velocity profile, be computed everywhere under the probe as in Section 4.3. Then, in a second step, this velocity profile must be used to solve the well-known advection–diffusion equation [12], which we shall admit here:


In typical MFP operations, the injected reagents in the gap between the surface and the mesa are typically of uniform concentration, that is, not varying in the z direction (). That observation is also strengthened as the gap G is most often kept small compared to other dimensions, making diffusion extremely fast in the vertical direction. A special case of MFP where this condition is broken is the hierarchical flow confinement MFP configuration [13]. Therefore, in most MFP cases, the advection–diffusion problem is also a quasi-2D problem arising from a quasi-2D flow profile. In such cases, the advection–diffusion equation takes a special form in which the velocity profile can be written as the gradient of the potential functions defined in Eqs. (4.19) and (4.20). Using the dimensionless groups used in Section 4.3, , , and , – and adding new groups related to diffusion, , , and , we obtain the normalized 2D advection–convection equation governing all MFP geometries:

Except for the simplest cases, mostly involving constant flows in shallow straight channels and the simplest boundary conditions, Eq. (4.26) is a challenging mathematical problem and often requires to be solved numerically using finite element softwares (see, e.g., Chapter 3). While numerical solutions are perfect to convey a general appreciation of the diffusion profile, they do not reveal the relationships (or scaling laws) between operation variables and observables. These relationships are crucial to provide MFP design and operation criteria. They are also essential for probe automation, as it is not possible to numerically compute the flow characteristics resulting from a certain set of operating conditions in real time during an experiment. A significant effort was dedicated to providing key scaling laws for diffusive processes in dipolar and quadrupolar MFP operations [10, 11, 14] and will be summarized here.

4.4.2 High Péclet Number Asymptotic Solutions Near Stagnation Points

Obtaining an exact solution of Eq. (4.26) in dipolar and quadrupolar MFP geometries is a challenging mathematical problem that yet remains to be solved. However, several approximations have been met to obtain effective models that are valid in some specific parts of the MFP domain and for specific Péclet numbers. They are extremely useful to determine the shape and “sharpness” of gradients at interfaces as well as to determine the extent of the blurring effect.

The approach used in this section consists of studying the diffusion locally in regions where the velocity is slowest, observing that these regions will possess the lowest Péclet numbers and thus the highest diffusion length scales. To do so, we linearize the velocity profiles generated from the potential functions in Eqs. (4.19) and (4.20) using a multivariable Taylor expansion of the second order. Then, we compute the local Péclet number and solve the linearized advection–diffusion locally. The values obtained will constitute upper bounds on the diffusion profiles. Values are then compared with full 3D FEM simulations to get a precise idea of the validity of the approximation. Not surprisingly, the position of the stagnation points in the flow patterns will play an important role here as these points are points of minimal velocity () and thus maximum diffusion. Floating Gradient Along the Central Line in a Microfluidic Quadrupole

When two different reagents (e.g., pure water () and an aqueous solution of small molecules ()) are injected by two opposing apertures in a microfluidic square quadrupole, the two flows impinge right under the center of the probe and make right angle turns before being re-aspirated by the aspiration apertures. This geometry creates a flat interface at the line of symmetry crossing the two aspiration apertures. Diffusion takes place right at the central stagnation point, and the diffused reagents are carried out rapidly to the aspiration apertures. This problem was first analyzed by Qasaimeh et al. [14] and, when taking the dominant terms in a multivariate Taylor expansion, yields a well-known flow configuration with an exact solution for high Péclet numbers. Such a Taylor expansion, keeping terms up to the second order in Eq. (4.20), evaluated near the central (inner) stagnation point () yields


where , in the mean flow rate of the quadrupole.

Using this result into Eq. (4.26) gives, in a dimensional form,

where . This problem has a general solution in terms of Hermite polynomials [15]. Yet, at high Péclet numbers and when considering the steady-state regime only, diffusion in the y direction can be neglected along the x-axis and the equation further reduces to


When Pe is sufficiently high, the diffusion length scale is much smaller than the distance between the injection apertures, and they can be construed as sources of constant concentration located infinitely far away, giving simple boundary conditions of the Dirichlet form . The solution to this 1D advection–diffusion problem for high Pe is of the form

The argument of the error function in Eq. (4.30) dictates the “sharpness” of the concentration gradient. The higher the Péclet number or the smaller the interaperture distance d, the sharper the resulting gradient. An experimentally measurable criterion for the gradient length can be defined by the distance over which reagent concentration goes from 10% to 90% of its maximum value. Numerically, this is given by

This criterion is compared with both numerical simulations and experimental results in Figure 4.6. For typical operating conditions using fluorescein, such as those from Qasaimeh et al. [14], with the values , , and , we get , and . Decreasing the above Pe value by a factor of 100 would increase Lgrad 10-fold to 370 µm. The high Pe solution regime is valid as long as . Experimental and numerical results from full FEM simulations demonstrate the accuracy of the analytical model even far away from the stagnation point. Thus, the simple result described in Eq. (4.31) can be used to generate precise concentration gradients by varying the mean flow rate and the gap, the only two operation parameters on which it depends. By increasing both Qinj and Qasp proportionally (keeping α constant) or by decreasing the gap G, it is thus possible to increase the gradient sharpness without altering the flow profile, that is, without altering the shape of the HFC, as the latter only depends on the interaperture distance d and the flow rate ratio α.

Figure 4.6 Comparison between advection–diffusion analytical models using the Hele-Shaw theory and full 3D simulations of advective–diffusive transport using FEM. Light gray: normalized gradient length as a function of flow rate ratio for microfluidic dipole probes with different aperture size ratios. Theory (solid line) compared with 16 different numerical results (markers). , , , , . Dark gray: normalized gradient lengths as a function of flow rate ratio in microfluidic quadrupole probes. Theory (solid line) compared with experimental data (markers). , , , , . Diffusion Broadening in the HFC Envelope for Dipoles and Quadrupoles

A similar approach can be used to compute the extent of the diffusion broadening of the HFC area outside dipoles and quadrupoles. This time, the velocity potential must be approximated in the vicinity of the outer stagnation points, whose positions are given by Eqs. (4.22) and (4.24). A multivariable Taylor expansion of the velocity potentials to the second order gives, for both dipoles and quadrupoles,


An analysis similar to the one presented in Section 4.2.1 holds here as well, using the same velocity profile functions yet with different expressions for the Péclet number (), . When making the approximation that the radius of curvature of the HFC is large compared with the extent of the diffusion broadening (valid again only at sufficiently high Pe), the analysis yields the same 1D steady-state differential equation (Eq. 4.28) with the same boundary conditions. In these cases, the solution is also of the form in Eq. (4.30) and the relative size of the diffusion broadening near the edges of the HFC is given by Eq. (4.31):

4.4.3 Numerical Investigation of Model Accuracy

In order to obtain a simple scaling law for the diffusion broadening related to a simple 1D model of the diffusion profile across the probe's HFC edges, a few approximations had to be made. Firstly, diffusion has to be a perturbation on the interfaces, such that the apertures can always be considered infinitely far away from the diffusive front (). Secondly, the interface radius of curvature is much larger than the extent of the diffusion length scale. Finally, the streamlines can be well approximated by a Hele-Shaw flow model.

These conditions are instantly met when the Péclet number is sufficiently high, providing a sharp, thin diffusion interface at the HFC edge, and the flow rate ratio is sufficiently low, reducing the interface curvature and preventing the stagnation points from being close or under the apertures. They are typically met in most applications of MFPs, such that the approximation is useful in practice. Figure 4.6 compares the 1D convection diffusion model described by Eqs. (4.30) and (4.34) with FEM solutions of the problem for the dipole and experimental results for the quadrupole.

For the particular case of the dipole, the large deviations observed between the simulation and the analysis at low α are due to the fact that, in previously published results, the low flow rates employed (), the Péclet number associated with is , thus breaking the high Pe condition. When α increases, Pe increases with and reaches , for . Yet the analytical model remains in agreement with the simulated data.

4.5 Conclusion

Since the advent of the concept of MFP in 2005, applications for this type of open microfluidic devices have multiplied. Their potential has now been demonstrated in various applications such as surface patterning and liquid-based microfabrication [9], controlled perfusion of brain tissue slices, gradient generation [16], selective cell extraction and processing in 2D cultures [17], and fixed tissue staining in immunohistochemistry [18]. The various probe concepts have also been diversified to meet the application needs, from simple push–pull configurations to quadrupole configurations performing gradient generation or complex nested reagent delivery and surface patterning, such as the hierarchical hydrodynamic flow confinement probe (hHFC) [13].

However, probe designs so far have always worked on the principle of injecting and aspirating fluids from apertures in a flat surface located in close proximity to a substrate to pattern or analyze. In this chapter, we provide a general approach, based on advection–diffusion Hele-Shaw flows, to model flow under MFPs of various shapes or forms. We have shown that the point source approximation for an aperture of any shape is very accurate, especially when apertures are relatively far away from each other. We have also shown how to analyze flow under an MFP from a distribution of point sources and how diffusive transport could be approximated using a simple 1D diffusion model developed near stagnation points (the points of maximum local Péclet numbers). Finally, we have discussed the strengths and limits of the Hele-Shaw model by comparing it with results from full 3D numerical simulations solving coupled Navier–Stokes and advection–diffusion equations. The main limitations of the Hele-Shaw model occur when the phenomena of interest occur at a scale that is on the same order or smaller than the thickness of the gap. But, as Henry Selby Hele-Shaw himself pointed it out in his seminal paper, regaining the Hele-Shaw features of a flow is always possible by decreasing the thickness of the cell (i.e., the gap between the probes in the surface in our case). In practice, this is not always possible or desirable, but the wide range of validity of the approach makes it widely applicable in any microfluidics problem when the flow domains are much shallower in one dimension compared with the two others.

Finally, the simple analytical or scaling models described in this chapter will find immediate application to automate probe operation, calibrate it, and systematize experimental design as suggested by Safavieh et al. [10] A recent example of the use of Hele-Shaw flows in probe operation and calibration is the floating MFP, where the gap height is fixed by a force balance model between probe weight and lift forces from the flow [19]. One can confidently expect that many other direct new MFP designs, using novel configurations and operation modes, will arise from an improved understanding of their behavior using the Hele-Shaw theory.


  1. 1 Hele-Shaw, H.S. (1898) The flow of water. Nature, 58, 34–36.
  2. 2 Bear, J. (1972) Dynamics of Fluids in Porous Media, Dover.
  3. 3 Gustafsson, B. and Vasil'ev, A. (2006) Conformal and Potential Analysis in Hele-Shaw cells, Birkhauser Verlag.
  4. 4 Kirby, B.J. (2010) Micro- and Nanoscale Fluid Mechanics: Transport in Microfluidic Devices, Cambridge University Press.
  5. 5 Batchelor, G.K. (2000) An Introduction to Fluid Dynamics, Cambridge University Press.
  6. 6 Bruus, H. (2008) Theoretical Microfluidics, Oxford University Press.
  7. 7 Christ, K.V. and Turner, K.T. (2011) Design of Hydrodynamically Confined Microfluidics: Controlling Flow Envelope and Pressure. Lab Chip, 11, 1491–1501.
  8. 8 Kaigala, G.V., Lovchik, R.D., Drechsler, U., and Delamarche, E. (2011) A vertical microfluidic probe. Langmuir, 27, 5686–5693.
  9. 9 Juncker, D., Schmid, H., and Delamarche, E. (2005) Multipurpose microfluidic probe. Nat. Mater., 4, 622–628.
  10. 10 Safavieh, M., Qasaimeh, M.A., Vakil, A., Juncker, D., and Gervais, T. (2015) Two-aperture microfluidic probes as flow dipole: theory and applications. Sci. Rep., 5, 11943.
  11. 11 Gervais, T., Safavieh, M., Qasaimeh, M.A., and Juncker, D. (2014) Systematic analysis of microfluidic probe design and operation. 36th Annual International Conference of the IEEE Engineering in Medicine and Biology, pp. 1567–1570.
  12. 12 Deen, W.M. (1998) Analysis of Transport Phenomena, Oxford.
  13. 13 Autebert, J., Kashyap, A., Lovchik, R.D., Delamarche, E., and Kaigala, G.V. (2014) Hierarchical hydrodynamic flow confinement: efficient use and retrieval of chemicals for microscale chemistry on surfaces. Langmuir, 30, 3640–3645.
  14. 14 Qasaimeh, M.A., Gervais, T., and Juncker, D. (2011) Microfluidic quadrupole and floating concentration gradient. Nat. Commun., 2, 464.
  15. 15 Lebedev, N.N. (1972) Special Functions and Their Applications, Dover.
  16. 16 Qasaimeh, M.A., Ricoult, S.G., and Juncker, D. (2013) Microfluidic probes for use in life sciences and medicine. Lab Chip, 13, 40–50.
  17. 17 Sarkar, A., Kolitz, S., Lauffenburger, D.A., and Han, J. (2014) Microfluidic probe for single-cell analysis in adherent tissue culture. Nat. Commun., 5, 3421.
  18. 18 Lovchik, R.D., Kaigala, G.V., Georgiadis, M., and Delamarche, E. (2012) Micro-immunohistochemistry using a microfluidic probe. Lab Chip, 12, 1040–1043.
  19. 19 Hitzbleck, M., Kaigala, G.V., Delamarche, E., and Lovchik, R.D. (2014) The floating microfluidic probe: distance control between probe and sample using hydrodynamic levitation. Appl. Phys. Lett., 104, 263501.