The deformation of microfluidic channels in a soft elastic medium has a central role in the operation of lab-on-a-chip devices, fluidic soft robots, liquid metal (LM) electronics, and other emerging soft-matter technologies. Understanding the influence of mechanical load on changes in channel cross section is essential for designing systems that either avoid channel collapse or exploit such collapse to control fluid flow and connectivity. In this paper, we examine the deformation of microchannel cross sections under far-field compressive stress and derive a “gauge factor” that relates externally applied pressure with change in cross-sectional area. We treat the surrounding elastomer as a Hookean solid and use two-dimensional plane strain elasticity, which has previously been shown to predict microchannel deformations that are in good agreement with experimental measurements. Numerical solutions to the governing Lamé (Navier) equations are found to match both the analytic solutions obtained from a complex stress function and closed-form algebraic approximations based on linear superposition. The application of this theory to soft microfluidics is demonstrated for several representative channel geometries.
Microfluidic systems for lab-on-a-chip diagnostics, biomaterials printing, and pathogen detection have had a transformative impact on medicine and biotechnology [1–4]. These devices are typically composed of microfluidic channels embedded in soft elastic materials and are produced using “soft lithography” fabrication techniques [5–7]. In recent years, soft microfluidic architectures have also been used for stretchable electronics with liquid metal (LM) [8–10], organ-on-a-chip tissue engineering [11,12], and fluid-powered soft robotics [13,14]. In most of these applications, the mechanical deformation of the embedded microfluidic channels and surrounding elastic medium are strongly coupled—leading to elasto-mechanical dependencies that can significantly influence device operation. Understanding the mechanics of these interactions is essential for developing predictive modeling tools that can inform the design and operation of soft microfluidic systems.
In many applications, elasto-mechanical coupling is an undesired consequence of using soft materials that can interfere with device functionality or performance. This is particularly true with stretchable electronics that use microfluidic LM wiring for soft, stretchable, and elastic circuits. For these circuits, microchannel deformation typically causes the electrical resistance to increase and can alter the passive electronic properties of the circuit . However, in some cases, elasto-mechanical coupling can be exploited to achieve unique functionalities that are not possible with a rigid device. This includes microfluidic “Quake” valves , LM pressure sensors and strain gauges [17–20], fluidic actuation [13,14], and pressure-controlled microfluidics that filter or capture nanoparticles [21,22] (Fig. 1(a)). In general, a theoretical model that relates the deformation of the microchannel and surrounding elastomer can inform device operation and enable the system to be designed in a way to either enhance or mitigate this coupling.
This paper presents theories that relate the influence of far-field stress on the deformation of a microchannel embedded within a soft elastic medium (Fig. 1(b)). Solutions are obtained for the cross-sectional geometries shown in Fig. 1(c), which are representative of the channel shapes typically used in applications. The elastic medium is treated as a Hookean solid with modulus E and Poisson's ratio ν. Previously it was shown that theoretical predictions based on Hooke's law are in strong agreement with experimental measurements for microchannels with rectangular  and triangular  cross section. The field equations and numerical solution method are presented in Sec. 2 along with exact and approximate analytic solutions based on classical solutions from the theory of linear elasticity. The numerical, exact, and approximate solutions are compared in Sec. 3 and shown to be in good agreement. This is followed by a discussion of how the proposed model can inform the design of selected systems. The paper closes with a brief summary of the results and prospects for future progress within the still nascent field of soft microfluidic mechanics and modeling (Sec. 4).
σ31 = σ13 = σ32 = σ23 = 0. Here, the notation applies to functions .
Here, corresponds to the top and bottom surfaces (X2 = ±L2/2) and t is the prescribed surface traction. In component form, on X2 = ±L2/2.
A numerical solution to the above boundary value problem is obtained in matlab R2015a (The Mathworks, Inc., Natick, MA) using the pdenonlin function. This solver uses the method of finite elements in which the domain is converted to a triangular mesh. The finite element matrices are then generated on the mesh and solved using a damped Gauss-Newton iteration algorithm. Numerical solutions are obtained for thin sheets with ν = 0.49 (i.e., incompressible) and unitless dimensions L1 = 10 and L2 = 1 and channels of width w = 0.1 and height h ranging from 0.01 to 0.1. If units of millimeters are selected, then the selected values are representative of typical microfluidic systems.
In soft microfluidics, the parameter is an important metric for elasto-mechanical coupling that relates the “sensitivity” of the microchannel geometry to an applied surface traction.
where , and are the complex conjugates of z, ϕ(z), and ψ(z), respectively. The prime in Eq. (12) denotes the derivative with respect to z.
Here, R is a scaling factor that depends upon the dimensions of the hole geometry. Since the series in Eq. (14) converges rapidly, we truncate the expression to three terms. For an ellipse of dimensions (a, b), , and for a rectangle of dimensions (w, h), , where and .
Finally, the displacement field is obtained by substituting Eqs. (14), (19), and (20) into Eq. (13). The area of the hole is numerically determined by importing the vector of contour points into the polyarea function in matlab 2016b. The gauge factor is calculated by . As discussed in Sec. 3, the values of A/A0 and are plotted in Figs. 2 and 3 and compared with the numerical and approximate analytic solutions.
The displacement field , microchannel area A, and gauge factor can be approximated by adapting solutions from classical problems in linear elasticity. As shown in Sec. 3, these closed-form approximations are found to be in reasonable agreement with the numerical and exact analytic solutions and can be useful for establishing design and operational guidelines for soft microfluidic systems.
Elliptical Cross Section.
We begin with the deformation of the elliptical cross section shown in Fig. 1(c). This can be determined by first examining the deformation of a circular cross section of radius R0 under pressure p0 such that the channel deforms into an ellipse with principle axes a and b (along the E1 and E2 directions, respectively). By linear superposition, the deformation of this ellipse under a surface traction p is equal to the deformation of a circular channel of radius R0 subject to an applied pressure of . Referring to Fig. 4, this corresponds to selecting different configurations for the natural () and reference placements () and decomposing the deformation mapping into an intermediate mapping and final mapping .
where α = h/w is the aspect ratio.
and the gauge factor is . As would be expected for a relatively small channel in a large linear elastic medium, resistance to deformation under a normalized pressure p/E decreases with an increase in the height-to-width aspect ratio of the cross section.
Rectangular Cross Section.
Again, the gauge factor decreases monotonically with increasing channel aspect ratio. As α gets large, converges to 1 − ν − 2ν2, which corresponds to the case of uniform strain (A ≈ Aε).
Results and Discussion
Numerical and analytic solutions for an elliptical cross section are presented in Fig. 2 for the following representative parameters: L1 = 10, L2 = 1, w = 0.2, h = 0.1, ν = 0.49. Figure 2(a) shows deformation under a normalized surface traction . The dashed lines correspond to the boundary of the unloaded state and the solid line corresponds to the numerical solution for the deformed shape obtained using the steps described in Sec. 2.1. The numerical solution for channel deformation appears to be in good agreement with the analytic approximations, shown as the dash-dot line. As shown in Fig. 2(b), there is also a reasonable agreement for calculations of the cross-sectional area as a function of applied pressure (). A comparison for the case of rectangular cross section is presented in Figs. 2(c) and 2(d).
As shown in Fig. 3, there also appears to be good agreement between numerical and analytic calculations of the gauge factor . For these plots, L1 = 10, L2 = 1, ν = 0.49, and the channel aspect ratio α ranges from 0.1 to 1. In the case of an elliptical cross section (Fig. 3(a)), the agreement appears to be close to exact. Figure 3(b) shows that Eq. (31) is in strong agreement with the numerical solution for a rectangular channel. As shown in Fig. 3(c), the gauge factor of microchannels with triangular (triangle markers) and diamond-shaped (diamond markers) cross section exhibit a similar monotonic dependency on aspect ratio. This can be accurately captured with the approximation , which appears to be in good agreement with the numerical solutions in the range of α ∈ [0.1, 1].
Both the numerical and analytic solutions are only valid when strain is small and the elastomer can be modeled as a linear elastic solid. For large strains, Hooke's law should be replaced with a nonlinear constitutive model (Neo-Hookean, Mooney-Rivlin, Ogden) and Cauchy stress should be balanced in the Eulerian description (i.e., ). For some cases, pdenonlin can still be used, but only with a tensor C = C(u, ∇u) that accounts for the stretch and rotation of the material coordinates. Otherwise, a more general finite element analysis (FEA) program is required.
where again, χ is the deformation mapping. The collapse of soft microfluidic channels under elastomer compression has previously been modeled using commercial FEA packages. This includes rectangular and triangular channels as well as triangular channels embedded with rigid spherical beads. While pressure-induced channel collapse is undesirable for LM circuits and sensors, it is central to the operation of soft microfluidic valves and particle filtration networks.
Hydrostatic pressure (ph) within the microfluidic channel can typically be ignored when surface pressure is concentrated over a small portion of the channel. For such cases, the reduction in internal volume from localized elastomer deformation is compensated by an equal increase distributed over a much larger volume (away from where pressure is applied). The latter requires only modest tractions on the channel wall and the corresponding hydrostatic pressure is expected to be small compared to the concentrated surface pressure. For cases when ph is significant, Eq. (3)2 should be replaced with the boundary condition . Likewise, analytic estimates may be obtained by superposing the approximations above with solutions for the case when traction is only applied on the channel walls.
The theoretical results and approximate models presented here have the potential to inform the choice of materials, designs, and operational conditions for a variety of soft microfluidic systems. In particular, we find that the gauge factor generally ranges from ∼1 to 10 for typical cross-sectional geometries. Therefore, in order for soft microfluidic systems to be resistant to collapse under applied pressure, the elastic modulus of the elastomer should be at least 100× greater than the maximum anticipated pressures. Alternatively, microfluidic channels that are designed to collapse under pressure should be embedded in an elastomer with a modulus that is approximately equal or 1/10th smaller in magnitude. This is especially useful in designing microfluidic channels for sensing, valving, or tunable particle filtration. For example, a circular channel (a = b) filled with ionically conductive fluid or liquid metal embedded in an incompressible elastomer has a gauge factor of 1.64. Based on Ohm's law, i.e., ΔR/R0 = (1 + ΔA/A0)−1 − 1, this implies that the electrical resistance R will increase by 50% if a pressure p = 0.2E is applied.
Finally, although only elastomer sheets with a single channel are modeled, the numerical approach can easily be extended to multiple embedded channels. For the case when channels are sparsely distributed (i.e., center-to-center spacing ≫ channel width), the solutions based on single channel deformation are expected to provide a reasonable approximation. This is because the influence of the channel on the internal strain field is highly localized and expected to have only modest influence on the mechanics of neighboring inclusions.
The compression of microchannels in an elastic medium has important implications in the design and operation of soft microfluidic systems. Of particular interest is the gauge factor that relates the relative change in channel area (ΔA/A0) with the normalized pressure (p/E) exerted on the surface of the surrounding elastomer. For the case of rectangular and elliptical cross section, there is good agreement between closed-form analytic approximations based on classical solutions in 2D elasticity with numerical solutions to the governing Lamé equations. In the case of a pressure-controlled microfluidic valve, filter, or sensor, the approximations can be useful for selecting material stiffness and channel aspect ratio for an anticipated range of pressures.
The numerical method presented here is applicable for examining small deformation of any prismatic channel geometry, including triangle, diamond, and other polygon-shaped cross section. For larger deformations, the theory must be modified to account for finite elastic deformation and unilateral contact between channel walls (i.e., collapse). Three-dimensional FEA modeling is required for accurate modeling of more general microfluidic networks with multiple channels and nonprismatic geometries. However, even for these general systems, the approximations presented here can help guide channel design and materials selection.
This work was supported by an Office of Naval Research BRC (Bio-Inspired Autonomous Systems, Dr. Tom McKenna).
Office of Naval Research (N00014-7-1-2063).
It should be noted that for problems in which large deformations are anticipated (finite elasticity), the elastic medium should be treated as a hyperelastic solid with an internal Cauchy stress that is divergence free in the spatial description, i.e., use the Eulerian nabla operator for ∇.
Also known as Navier's or Lamé's equations.
Otherwise, the coefficient 1/2 in Eq. (9) should be replaced with –1/2.