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.

## Introduction

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 [15]. 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 [16], 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 [18] and triangular [20] 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).

## Theory

*L*

_{1},

*L*

_{2},

*L*

_{3}} is embedded with a single microchannel of width

*w*and height

*h*that extends along its entire depth. Referring to the coordinates in Fig. 1(b), the sheet is subject to compressive pressure

*p*on its top (

*X*

_{2}=

*L*

_{2}/2) and bottom (

*X*

_{2}= –

*L*

_{2}/2) surfaces. This induces a displacement field

**u**=

*u*(

_{i}**X**)

**E**

*and stress $\sigma =\sigma ij(X)Ei\u2297Ej$ within the elastomer. Here, $X=XiEi\u2208B0,\u2009B0$ is the natural placement of the elastic medium in Euclidean space ($E$), and the indices*

_{i}*i*,

*j*∈ {1, 2, 3} are subject to Einstein summation convention. Assuming plane strain loading, $\sigma $ has components [24]

*σ*_{31} = *σ*_{13} = *σ*_{32} = *σ*_{23} = 0. Here, the notation $f,i=\u2202f/\u2202Xi$ applies to functions $f:B0\u2192\mathbb{R}$.

^{3}Letting $\psi =(1\u22122\nu )/2(1\u2212\nu )$ and

*ν*

_{1}=

*ν*/(1 −

*ν*), it follows from Eq. (1) that the balance law can be expressed by the following pair of partial differential equations:

Here, $B0t$ corresponds to the top and bottom surfaces (*X*_{2} = ±*L*_{2}/2) and **t** is the prescribed surface traction. In component form, $t=\u2213pE2$ on *X*_{2} = ±*L*_{2}/2.

### Numerical Solution.

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 $B0$ 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 *L*_{1} = 10 and *L*_{2} = 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.

*pdenonlin*as a system of elliptic partial differential equations with the general form

**C**,

**a**, and

**f**are fourth-, second-, and first-order tensors, respectively. For the current problem, $\u2207=(\u2202/\u2202X1)E1+(\u2202/\u2202X2)E2,\u2009a=f=0$, and $C=CijklEi\u2297Ej\u2297Ek\u2297El$. From Eq. (2), it follows that

**C**are zero. In matlab,

**C**is imported as a column vector of character strings with elements (1, 0, 0,

*ψ*, 0,

*ψ*,

*ν*

_{1}, 0, 0,

*ν*

_{1},

*ψ*, 0,

*ψ*, 0, 0, 1). Likewise, the boundary conditions in Eq. (3) must be expressed as

*X*

_{2}= ±

*L*

_{2}/2. For convenience, we will also define a normalized surface pressure $p\u0302=p/E$. As shown next, the normalized surface pressure $p\u0302$ will also be useful in calculating the

*gauge factor*that relates applied pressure with the change in cross-sectional area Δ

*A*of the microchannel.

*x*,

^{c}*y*) are obtained from the deformation mapping

^{c}**x**=

**X**+

**u**for points $Xc\u2208C$. In parametric form, $X\alpha c=X\alpha c(t)$ for

*α*∈ {1, 2}, where

*t*can be an arclength along the boundary or angle defined with respect to the

*X*

_{1}axis. The parameter

*t*must increase when moving along $C$ in the counter-clockwise direction.

^{5}In discrete form where

*t*is an array of length

*n*, the contour is discretized by points $xic=(X1)ic+u1(XiC)$ and $yic=(X2)ic+u2(Xic)$, where $(X\alpha )ic=X\alpha c(ti)$ and

*i*∈ {1,…,

*n*}. The area is calculated as

In soft microfluidics, the parameter $G$ is an important metric for elasto-mechanical coupling that relates the “sensitivity” of the microchannel geometry to an applied surface traction.

### Analytic Solution.

*x*

_{1},

*x*

_{2}) within the elastic media by a complex variable,

*z*=

*x*

_{1}+

*ix*

_{2}. Using the Kolosov-Muskhelishvili formulae [24,25], the stress and displacement fields can be expressed in terms of analytic complex stress potential functions

*ϕ*(

*z*),

*ψ*(

*z*)

where $\kappa =3\u22124\nu ;\u2009z\xaf,\u2009\varphi (z)\xaf$, and $\psi (z)\xaf$ are the complex conjugates of *z*, *ϕ*(*z*), and *ψ*(*z*), respectively. The prime in Eq. (12) denotes the derivative with respect to *z*.

*z*=

*ω*(

*ζ*), is employed to map the region in $B0$ outside the hole, from the

*z*plane to the interior of a unit circle in the

*ζ*plane

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*), $m1=(a\u2212b)/(a+b),\u2009m2=0,\u2009m3=0$, and for a rectangle of dimensions (*w*, *h*), $m1=(a+a\xaf)/2,\u2009m2=(a\u2212a\xaf)2/24,\u2009m3={(a2\u2212a\xaf2)(a\u2212a\xaf)}/80$, where $a=e2ik\pi $ and $k=0.068\u2009log(38.314\xb7h/w)$.

*ϕ*(

*ζ*) and

*ψ*(

*ζ*) are found by applying Cauchy integrals to the derived boundary conditions on the contour

*γ*of the fictitious hole in the parametric

*ζ*-plane and solving the resulting pair of equations [27]

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 $G=(1\u2212A/A0)/p\u0302$. As discussed in Sec. 3, the values of *A*/*A*_{0} and $G$ are plotted in Figs. 2 and 3 and compared with the numerical and approximate analytic solutions.

### Approximate Solution.

The displacement field $u:B0\u2192B$, microchannel area *A*, and gauge factor $G$ 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 *R*_{0} under pressure *p*_{0} such that the channel deforms into an ellipse with principle axes *a* and *b* (along the **E**_{1} and **E**_{2} 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 *R*_{0} subject to an applied pressure of $p\u2032=p0+p$. Referring to Fig. 4, this corresponds to selecting different configurations for the natural ($B0\u2032$) and reference placements ($B0$) and decomposing the deformation mapping $\chi \u2032=\chi \u2218\chi 0$ into an intermediate mapping $\chi 0:B0\u2032\u2192B0$ and final mapping $\chi :B0\u2192B$.

*Michell problem*[24,25]. For a compressive stress

*p*

_{0}, the displacement field around the hole boundary has components

*φ*= (

*k*+ 1)/8,

*k*= 3 − 4

*ν*, and the polar coordinate

*θ*is defined with respect to the

*X*

_{1}axis. To achieve an elliptic cross section with dimensions

*a*=

*w*/2 and

*b*=

*h*/2, the elastic medium must be deformed such that

*u*

_{1}(0) =

*a*−

*R*

_{0}and

*u*

_{2}(

*π*/2) =

*b*−

*R*

_{0}. This implies the following opening radius and surface pressure:

where *α* = *h*/*w* is the aspect ratio.

and the gauge factor is $Ge=(1\u2212\nu 2)(9\u2212\alpha 2)/4\alpha $. 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.

*A*/

*A*

_{0}and $G$ can be obtained by examining the dual case of tensile loading. Under tension, the final area can be estimated by superposing the area $A\epsilon =wh(1+\epsilon )(1\u2212\nu (1+\nu )\epsilon )$ under a uniform global (average) strain

*ε*=

*σ*(1 −

*ν*

^{2})/

*E*with the increase in area (Δ

*A*) for a

*crack opening*under a far-field tension

*σ*. The latter is obtained from linear elastic fracture mechanics, by treating the channel cross section as a narrow slit of width

*w*. Referring to Fig. 4, the vertical displacement of the slit boundary (

*u*

_{2}) is obtained for a stress intensity factor $KI=\sigma \pi w/2$ [24,29,30]

*σ*is replaced with –

*p*and the position of points along the boundary is estimated as

Again, the gauge factor decreases monotonically with increasing channel aspect ratio. As *α* gets large, $Gp$ 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: *L*_{1} = 10, *L*_{2} = 1, *w* = 0.2, *h* = 0.1, *ν* = 0.49. Figure 2(a) shows deformation under a normalized surface traction $|t\u0302|=(1\u2212\nu 2)p\u0302=0.05$. 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 ($0\u2264p\u0302\u22640.1$). 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 $G$. For these plots, *L*_{1} = 10, *L*_{2} = 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 $G=1.2/\alpha $, 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., $\u2207E\xb7\sigma =0$). 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.

*E*. This is especially true for channels with a low cross-sectional aspect ratio (

*α*). For such deformations, analysis should also include a

*unilateral constraint*to prevent interpenetration

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 (*p _{h}*) 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

*p*is significant, Eq. (3)

_{h}_{2}should be replaced with the boundary condition $\sigma \xb7n=\u2212phn\u2009\u2200X\u2208\u2202B0\\u2202B0t$. 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*/*R*_{0} = (1 + Δ*A*/*A*_{0})^{−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.

## Conclusion

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 $G$ that relates the relative change in channel area (Δ*A*/*A*_{0}) 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.

## Acknowledgment

This work was supported by an Office of Naval Research BRC (Bio-Inspired Autonomous Systems, Dr. Tom McKenna).

## Funding Data

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.