Deformation of Microchannels Embedded in an Elastic Medium

Ramachandran, Vivek; Majidi, Carmel

doi:10.1115/1.4040477

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.

Fig. 1

View large Download slide

(a) Examples of soft microfluidics for emerging applications: (i) Quake valve (reproduced with permission from Unger et al. [16] Copyright 2000 by AAAS), (ii) fluidic soft robot (reproduced with permission from Wehner et al. [14] Copyright 2016 by Springer Nature), (iii) LM electronics (reproduced with permission from Ozutemiz et al. [23] Copyright 2018 by Wiley—VCH Verlag GmbH & Co. KGaA), and (iv) tunable nanoparticle filtration (reproduced with permission from Huh et al. [21] Copyright 2007 by Springer Nature). (b) Soft microfluidic channel modeled as a prismatic opening embedded inside of an elastic sheet subject to surface pressure p. (c) Selected channel cross section, which is representative of geometries typically used in soft microfluidics.

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

A rectangular elastic sheet of dimensions {L₁, L₂, L₃} 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₂ = L₂/2) and bottom (X₂ = –L₂/2) surfaces. This induces a displacement field u = u_i(X)E_i and stress

σ = σ_{i j} (X) E_{i} \otimes E_{j}

within the elastomer. Here,

X = X_{i} E_{i} \in B_{0}, B_{0}

is the natural placement of the elastic medium in Euclidean space (⁠

E

⁠), and the indices i, j ∈ {1, 2, 3} are subject to Einstein summation convention. Assuming plane strain loading,

σ

has components [24]

\begin{matrix} σ_{11} = \frac{E {(1 - ν) u_{1, 1} + ν u_{2, 2}}}{(1 + ν) (1 - 2 ν)} σ_{12} = σ_{21} = \frac{E (u_{1, 2} + u_{2, 1})}{2 (1 + ν)} \\ σ_{22} = \frac{E {ν u_{1, 1} + (1 - ν) u_{2, 2}}}{(1 + ν) (1 - 2 ν)} σ_{33} = \frac{E ν (u_{1, 1} + u_{2, 2})}{(1 + ν) (1 - 2 ν)} \end{matrix}

(1)

σ₃₁ = σ₁₃ = σ₃₂ = σ₂₃ = 0. Here, the notation $f_{, i} = \partial f / \partial X_{i}$ applies to functions $f : B_{0} \to ℝ$ ⁠.

For a linear elastic solid,

\nabla \cdot σ = 0

at static equilibrium, where

\nabla = {\partial / \partial X_{i}} E_{i}

is the Lagrangian (material) nabla operator.³ Letting

ψ = (1 - 2 ν) / 2 (1 - ν)

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

\begin{matrix} u_{1, 11} + ν_{1} u_{2, 21} + ψ (u_{1, 22} + u_{2, 12}) = 0 \\ u_{2, 22} + ν_{1} u_{1, 12} + ψ (u_{1, 21} + u_{2, 11}) = 0 \end{matrix}

(2)

The solution to this displacement formulation⁴ must satisfy the following boundary conditions:

σ \cdot n = t \forall X \in \partial B_{0}^{t} σ \cdot n = 0 \forall X \in \partial B_{0} \ \partial B_{0}^{t}

(3)

Here, $B_{0}^{t}$ corresponds to the top and bottom surfaces (X₂ = ±L₂/2) and t is the prescribed surface traction. In component form, $t = \mp p E_{2}$ on X₂ = ±L₂/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 $B_{0}$ 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₁ = 10 and L₂ = 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.

The governing equations in are imported into pdenonlin as a system of elliptic partial differential equations with the general form

\nabla \cdot (C : \nabla u) = a u + f

(4)

where C, a, and f are fourth-, second-, and first-order tensors, respectively. For the current problem,

\nabla = (\partial / \partial X_{1}) E_{1} + (\partial / \partial X_{2}) E_{2}, a = f = 0

⁠, and

C = C_{i j k l} E_{i} \otimes E_{j} \otimes E_{k} \otimes E_{l}

⁠. From Eq. , it follows that

\begin{matrix} C_{1111} = C_{2222} = 1 C_{1122} = C_{2211} = ν_{1} \\ C_{1212} = C_{1221} = C_{2112} = C_{2121} = ψ \end{matrix}

(5)

and the remaining elements of C are zero. In matlab, C is imported as a column vector of character strings with elements (1, 0, 0, ψ, 0, ψ, ν₁, 0, 0, ν₁, ψ, 0, ψ, 0, 0, 1). Likewise, the boundary conditions in Eq. must be expressed as

(C : \nabla u) \cdot n + q u = g

(6)

where q is a second-order tensor and g is the prescribed surface traction. Here, q = 0 and

g = {\begin{matrix} \hat{t} & \forall X \in \partial B_{0}^{t} \\ 0 & \forall X \in \partial B_{0} \ \partial B_{0}^{t} \end{matrix}

(7)

where the nondimensional traction can be expressed as

\hat{t} = \mp \frac{(1 + ν) (1 - 2 ν) p}{(1 - ν) E} E_{2}

(8)

for X₂ = ±L₂/2. For convenience, we will also define a normalized surface pressure

\hat{p} = p / E

⁠. As shown next, the normalized surface pressure

\hat{p}

will also be useful in calculating the gauge factor that relates applied pressure with the change in cross-sectional area ΔA of the microchannel.

The cross-sectional area of the microchannel is determined by using Green's theorem

A = \frac{1}{2} \oint_{C} {x^{c} d y^{c} - y^{c} d x^{c}}

(9)

Here, the contour

C

represents the space curve formed by the microchannel boundary in the natural placement (undeformed configuration). The coordinates (x^c, y^c) are obtained from the deformation mapping x = X + u for points

X^{c} \in C

⁠. In parametric form,

X_{α}^{c} = X_{α}^{c} (t)

for α ∈ {1, 2}, where t can be an arclength along the boundary or angle defined with respect to the X₁ axis. The parameter t must increase when moving along

C

in the counter-clockwise direction.⁵ In discrete form where t is an array of length n, the contour is discretized by points

x_{i}^{c} = {(X_{1})}_{i}^{c} + u_{1} (X_{i}^{C})

and

y_{i}^{c} = {(X_{2})}_{i}^{c} + u_{2} (X_{i}^{c})

⁠, where

{(X_{α})}_{i}^{c} = X_{α}^{c} (t_{i})

and i ∈ {1,…, n}. The area is calculated as

A = \frac{1}{2} \sum_{i = 2}^{n} {x_{i}^{c} Δ y_{i}^{c} - y_{i}^{c} Δ x_{i}^{c}}

(10)

where

Δ x_{i}^{c} = x_{i}^{c} - x_{i - 1}^{c}

and

Δ y_{i}^{c} = y_{i}^{c} - y_{i - 1}^{c}

⁠. Finally, the gauge factor corresponds to the ratio of the relative change in cross-sectional area ΔA/A₀ and normalized pressure

G = \lim_{p \to 0} \frac{1}{\hat{p}} {1 - \frac{A}{A_{0}}}

(11)

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.

As with other classes of problems in two-dimensional (2D) linear elasticity, an analytic solution to the above boundary value problem can be obtained using complex analysis. This involves replacing the coordinates (x₁, x₂) within the elastic media by a complex variable, z = x₁ + ix₂. 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)

\begin{matrix} σ_{11} + σ_{22} = 2 [ϕ'_{1} (z) + \bar{ϕ'_{1} (z)}] \\ σ_{22} - σ_{11} + 2 i σ_{12} = 2 [\bar{z} ϕ''_{1} (z) + ψ'_{1} (z)] \\ 2 μ (u_{1} + i u_{2}) = κ ϕ (z) - z \bar{ϕ' (z)} - \bar{ψ (z)} \end{matrix}

(12)

where $κ = 3 - 4 ν; \bar{z}, \bar{ϕ (z)}$ ⁠, and $\bar{ψ (z)}$ are the complex conjugates of z, ϕ(z), and ψ(z), respectively. The prime in Eq. (12) denotes the derivative with respect to z.

To find the displacement field around holes in an infinite media, a conformal transformation z = ω(ζ), is employed to map the region in

B_{0}

outside the hole, from the z plane to the interior of a unit circle in the ζ plane

2 μ (u_{1} + i u_{2}) = κ ϕ (ζ) - \frac{ω (ζ)}{\bar{ω' (ζ)}} \bar{ϕ' (ζ)} - \bar{ψ (ζ)}

(13)

where

ψ (ζ) = ψ \circ ω (ζ)

and

ϕ (ζ) = ϕ \circ ω (ζ)

⁠. The mapping function

ω (ζ)

takes the general form [26]

ω (ζ) = R [\frac{1}{ζ} + \sum_{k = 1}^{\infty} m_{k} ζ^{2 k - 1}]

(14)

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), $m_{1} = (a - b) / (a + b), m_{2} = 0, m_{3} = 0$ ⁠, and for a rectangle of dimensions (w, h), $m_{1} = (a + \bar{a}) / 2, m_{2} = {(a - \bar{a})}^{2} / 24, m_{3} = {(a^{2} - {\bar{a}}^{2}) (a - \bar{a})} / 80$ ⁠, where $a = e^{2 i k π}$ and $k = 0.068 log (38.314 \cdot h / w)$ ⁠.

As reported in Ref. [25], the stress potential functions

{ϕ (ζ), ψ (ζ)}

can each be resolved into the pair of functions

\begin{matrix} ϕ_{1} (ζ) = \frac{p ω (ζ)}{4} \\ ψ_{1} (ζ) = \frac{p ω (ζ)}{2} \end{matrix}

(15)

and the following pair of holomorphic functions in the region

| ζ | < 1

⁠:

\begin{matrix} ϕ_{0} (ζ) = \sum_{n = 1}^{\infty} a_{n} ζ^{n} \\ ψ_{0} (ζ) = \sum_{n = 0}^{\infty} b_{n} ζ^{n} \end{matrix}

(16)

The expressions for ϕ(ζ) 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]

\begin{matrix} ϕ_{0} (ζ) + \frac{1}{2 π i} \int_{γ} \frac{ω (σ)}{\bar{ω' (σ)}} \frac{\bar{ϕ'_{0} (σ)}}{σ - ζ} d σ = \frac{1}{2 π i} \int_{γ} \frac{f_{1}^{0} + i f_{2}^{0}}{σ - ζ} d σ \\ ψ_{0} (ζ) + \frac{1}{2 π i} \int_{γ} \frac{\bar{ω (σ)}}{ω' (σ)} \frac{ϕ'_{0} (σ)}{σ - ζ} d σ = \frac{1}{2 π i} \int_{γ} \frac{f_{1}^{0} - i f_{2}^{0}}{σ - ζ} d σ \end{matrix}

(17)

The derived boundary conditions for the basic stress state of this medium are [25]

\begin{matrix} f_{1}^{0} + i f_{2}^{0} = \frac{p R}{2} [ω (ζ) + \bar{ω (ζ)}] \\ f_{1}^{0} - i f_{2}^{0} = \frac{p R}{2} [\bar{ω (ζ)} + ω (ζ)] \end{matrix}

(18)

As an example, the stress potential functions for an ellipse

{ϕ_{e} (ζ), ψ_{e} (ζ)}

and a rectangle

{ϕ_{r} (ζ), ψ_{r} (ζ)}

of aspect ratio

h / w = b / a = 0.5

are given by

\begin{matrix} ϕ_{e} (ζ) = \hat{p} R (\frac{1}{4 ζ} - 0.583 ζ) \\ ψ_{e} (ζ) = \hat{p} R (\frac{1}{2 ζ} + 1.5 ζ + \frac{6.67 ζ}{ζ^{2} - 3}) \end{matrix}

(19)

and

\begin{matrix} ϕ_{r} (ζ) = \hat{p} R (\frac{1}{4 ζ} - 0.486 ζ + 0.0382 ζ^{3}) \\ ψ_{r} (ζ) = \hat{p} R (\frac{1}{2 ζ} - \frac{0.773 ζ^{3} + 1.332 ζ}{0.918 ζ^{4} - 0.574 ζ^{2} + 2}) \end{matrix}

(20)

respectively. The equations in match those provided in previous studies [28].

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 - A / A_{0}) / \hat{p}$ ⁠. As discussed in Sec. 3, the values of A/A₀ and $G$ are plotted in Figs. 2 and 3 and compared with the numerical and approximate analytic solutions.

Fig. 2

View large Download slide

(a) Deformation of an elliptical channel with w = 0.2, h = 0.1, p = 0.05, and ν = 0.49: (gray) initial shape, (markers) numerical solution for deformed shape, (dashed) approximate solution for deformed shape. (b) Cross-sectional area as a function of normalized pressure: (circular markers) numerical solution, (dashed) approximate solution, and (solid) analytic solution. (c) Deformation of a rectangular channel with w = 0.2, h = 0.1, p = 0.05, and ν = 0.49: (gray) initial shape, (markers) numerical solution for deformed shape, (dashed) approximate solution for deformed shape. (d) Cross-sectional area as a function of normalized pressure: (circular markers) numerical solution, (dashed) approximate solution, and (solid) analytic solution.

Fig. 3

View large Download slide

Gauge factor $G$ versus channel aspect ratio α. (a) Circular cross section with a = 0.1 and varying b: (circular markers) numerical solution, (solid) analytic solution, (dashed) approximate solution. (b) Rectangular cross section with w = 0.1 and varying h: (square markers) numerical solution, (solid) analytic solution, and (dashed) approximate solution. (c) Numerical solutions for triangular (triangle markers) and diamond-shape (diamond markers) also exhibit a similar monotonic dependency; the dashed line corresponds to a fit of $G = 1.2 / α$ ⁠.

Approximate Solution.

The displacement field $u : B_{0} \to B$ ⁠, 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₀ under pressure p₀ such that the channel deforms into an ellipse with principle axes a and b (along the E₁ and E₂ 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₀ subject to an applied pressure of $p^{'} = p_{0} + p$ ⁠. Referring to Fig. 4, this corresponds to selecting different configurations for the natural (⁠ $B_{0}^{'}$ ⁠) and reference placements (⁠ $B_{0}$ ⁠) and decomposing the deformation mapping $χ^{'} = χ \circ χ_{0}$ into an intermediate mapping $χ_{0} : B_{0}^{'} \to B_{0}$ and final mapping $χ : B_{0} \to B$ ⁠.

Fig. 4

View large Download slide

(a) The solution to the classical Michell problem of a circular opening in an elastic plate under far-field stress is used to estimate the deformation of an elliptical channel with principle dimensions a and b. (b) The deformation of a rectangular channel is estimated by superposing uniform deformation of a homogenous Hookean solid with the opening of a slit under far-field stress.

Study of the deformation of a circular hole in an elastic medium with far-field uniaxial stress is classically known as the Michell problem [24,25]. For a compressive stress p₀, the displacement field around the hole boundary has components

u_{1} (θ) = \frac{2 (1 + ν)}{E} φ p_{0} R_{0} cos θ

(21)

u_{2} (θ) = - \frac{6 (1 + ν)}{E} φ p_{0} R_{0} sin θ

(22)

where φ = (k + 1)/8, k = 3 − 4ν, and the polar coordinate θ is defined with respect to the X₁ 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₁(0) = a − R₀ and u₂(π/2) = b − R₀. This implies the following opening radius and surface pressure:

p_{0} = \frac{E (1 - α)}{2 φ (1 + ν) (3 + α)}

(23)

R_{0} = \frac{w}{2} {1 + \frac{2 (1 + ν) φ p_{0}}{E}}^{- 1}

(24)

where α = h/w is the aspect ratio.

The final shape of the channel is determined by superposition

x_{1}^{c} = R_{0} cos θ {1 + \frac{2 φ (1 + ν) (p + p_{0})}{E}}

(25)

x_{2}^{c} = R_{0} sin θ {1 - \frac{6 φ (1 + ν) (p + p_{0})}{E}}

(26)

The corresponding cross-sectional area is

A = π R_{0}^{2} {1 + \frac{2 φ (1 + ν) (p + p_{0})}{E}} {1 - \frac{6 φ (1 + ν) (p + p_{0})}{E}}

(27)

and the gauge factor is $G_{e} = (1 - ν^{2}) (9 - α^{2}) / 4 α$ ⁠. 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.

For a rectangular cross section, closed-form algebraic approximations for A/A₀ 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_{ε} = w h (1 + ε) (1 - ν (1 + ν) ε)

under a uniform global (average) strain ε = σ(1 − ν²)/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₂) is obtained for a stress intensity factor

K_{I} = σ \sqrt{π w / 2}

[24,29,30]

u_{2} = \frac{4 (1 - ν^{2}) K_{I}}{E} \sqrt{\frac{w - 2 X_{1}}{4 π}}

(28)

Integrating over the length of the slit yields

Δ A = 4 (1 - ν^{2}) \sqrt{2} σ w^{2} / 3 E

⁠. For compression, σ is replaced with –p and the position of points along the boundary is estimated as

\begin{matrix} x^{c} = (1 + ν (1 + ν) \frac{p}{E}) X_{1}^{c} E_{1} + {(1 - (1 - ν^{2}) \frac{p}{E}) X_{2}^{c} \\ - sgn (X_{2}) \frac{4 (1 - ν^{2}) p \sqrt{w}}{2 E} \sqrt{\frac{w}{2} - X_{1}^{c}}} E_{2} \end{matrix}

(29)

The final cross-sectional area A = A_ε + ΔA is calculated as

\begin{matrix} A = w h (1 - (1 - ν^{2}) \frac{p}{E}) (1 + ν (1 + ν) \frac{p}{E}) \\ - \frac{4 (1 - ν^{2}) \sqrt{2} p w^{2}}{3 E} \end{matrix}

(30)

and the gauge factor is

G_{p} = 1 - ν - 2 ν^{2} + \frac{4 (1 - ν^{2}) \sqrt{2}}{3 α}

(31)

Again, the gauge factor decreases monotonically with increasing channel aspect ratio. As α gets large, $G_{p}$ converges to 1 − ν − 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₁ = 10, L₂ = 1, w = 0.2, h = 0.1, ν = 0.49. Figure 2(a) shows deformation under a normalized surface traction $| \hat{t} | = (1 - ν^{2}) \hat{p} = 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 \leq \hat{p} \leq 0.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₁ = 10, L₂ = 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 / α$ ⁠, 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., $\nabla_{E} \cdot σ = 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.

In the current deformation analysis, it is also assumed that the walls of the channel do not make contact or interpenetrate. In general, the channel will collapse as the surface pressure approaches 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

\exists! X_{1}, X_{2} \in \partial B_{0} \ C : χ (X_{1}) \equiv χ (X_{2})

(32)

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_h is significant, Eq. (3)₂ should be replaced with the boundary condition $σ \cdot n = - p_{h} n \forall X \in \partial B_{0} \ \partial B_{0}^{t}$ ⁠. 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₀ = (1 + ΔA/A₀)⁻¹ − 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₀) 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).

3

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

4

Also known as Navier's or Lamé's equations.

5

Otherwise, the coefficient 1/2 in Eq. (9) should be replaced with –1/2.

References

1.

Stone

,

H. A.

,

Stroock

,

A. D.

, and

Ajdari

,

A.

,

2004

, “

Engineering Flows in Small Devices: Microfluidics Toward a Lab-on-a-Chip

,”

Annu. Rev. Fluid Mech.

,

36

(

1

), pp.

381

–

411

.

Google Scholar

Crossref

2.

Squires

,

T. M.

, and

Quake

,

S. R.

,

2005

, “

Microfluidics: Fluid Physics at the Nanoliter Scale

,”

Rev. Mod. Phys.

,

77

(

3

), p.

977

.

Google Scholar

Crossref

3.

Whitesides

,

G. M.

,

2006

, “

The Origins and the Future of Microfluidics

,”

Nature

,

442

(

7101

), pp.

368

–

373

.

Google Scholar

Crossref

PubMed

4.

Sackmann

,

E. K.

,

Fulton

,

A. L.

, and

Beebe

,

D. J.

,

2014

, “

The Present and Future Role of Microfluidics in Biomedical Research

,”

Nature

,

507

(

7491

), pp.

181

–

189

.

Google Scholar

Crossref

PubMed

5.

Xia

,

Y.

, and

Whitesides

,

G. M.

,

1998

, “

Soft Lithography

,”

Annu. Rev. Mater. Sci.

,

28

(

1

), pp.

153

–

184

.

Google Scholar

Crossref

6.

Quake

,

S. R.

, and

Scherer

,

A.

,

2000

, “

From Micro-to Nanofabrication With Soft Materials

,”

Science

,

290

(

5496

), pp.

1536

–

1540

.

Google Scholar

Crossref

PubMed

7.

Qin

,

D.

,

Xia

,

Y.

, and

Whitesides

,

G. M.

,

2010

, “

Soft Lithography for Micro-and Nanoscale Patterning

,”

Nat. Protoc.

,

5

(

3

), pp.

491

–

502

.

Google Scholar

Crossref

PubMed

8.

Dickey

,

M. D.

,

Chiechi

,

R. C.

,

Larsen

,

R. J.

,

Weiss

,

E. A.

,

Weitz

,

D. A.

, and

Whitesides

,

G. M.

,

2008

, “

Eutectic Gallium-Indium (Egain): A Liquid Metal Alloy for the Formation of Stable Structures in Microchannels at Room Temperature

,”

Adv. Funct. Mater.

,

18

(

7

), pp.

1097

–

1104

.

Google Scholar

Crossref

9.

Cheng

,

S.

, and

Wu

,

Z.

,

2012

, “

Microfluidic Electronics

,”

Lab Chip

,

12

(

16

), pp.

2782

–

2791

.

Google Scholar

Crossref

PubMed

10.

Joshipura

,

I. D.

,

Ayers

,

H. R.

,

Majidi

,

C.

, and

Dickey

,

M. D.

,

2015

, “

Methods to Pattern Liquid Metals

,”

J. Mater. Chem. C

,

3

(

16

), pp.

3834

–

3841

.

Google Scholar

Crossref

11.

Huh

,

D.

,

Kim

,

H. J.

,

Fraser

,

J. P.

,

Shea

,

D. E.

,

Khan

,

M.

,

Bahinski

,

A.

,

Hamilton

,

G. A.

, and

Ingber

,

D. E.

,

2013

, “

Microfabrication of Human Organs-on-Chips

,”

Nat. Protoc.

,

8

(

11

), pp.

2135

–

2157

.

Google Scholar

Crossref

PubMed

12.

Bhatia

,

S. N.

, and

Ingber

,

D. E.

,

2014

, “

Microfluidic Organs-on-Chips

,”

Nat. Biotechnol.

,

32

, pp. 760–772.

13.

Wakimoto

,

S.

,

Ogura

,

K.

,

Suzumori

,

K.

, and

Nishioka

,

Y.

,

2009

, “

Miniature Soft Hand With Curling Rubber Pneumatic Actuators

,”

IEEE International Conference on Robotics and Automation

(

ICRA'09

), Kobe, Japan, May 12–17, pp.

556

–

561

.

14.

Wehner

,

M.

,

Truby

,

R. L.

,

Fitzgerald

,

D. J.

,

Mosadegh

,

B.

,

Whitesides

,

G. M.

,

Lewis

,

J. A.

, and

Wood

,

R. J.

,

2016

, “

An Integrated Design and Fabrication Strategy for Entirely Soft, Autonomous Robots

,”

Nature

,

536

(

7617

), pp.

451

–

455

.

Google Scholar

Crossref

PubMed

15.

Kim

,

H.-J.

,

Son

,

C.

, and

Ziaie

,

B.

,

2008

, “

A Multiaxial Stretchable Interconnect Using Liquid-Alloy-Filled Elastomeric Microchannels

,”

Appl. Phys. Lett.

,

92

(

1

), p.

011904

.

Google Scholar

Crossref

16.

Unger

,

M. A.

,

Chou

,

H.-P.

,

Thorsen

,

T.

,

Scherer

,

A.

, and

Quake

,

S. R.

,

2000

, “

Monolithic Microfabricated Valves and Pumps by Multilayer Soft Lithography

,”

Science

,

288

(

5463

), pp.

113

–

116

.

Google Scholar

Crossref

PubMed

17.

Whitney

,

R.

,

1949

, “

The Measurement of Changes in Human Limb-Volume by Means of a Mercury-Inrubber Strain Gauge

,”

J. Physiol.

,

109

(

1–2

), p.

Proc–5

.

18.

Majidi

,

C.

,

Park

,

Y.-L.

,

Kramer

,

R.

,

Bérard

,

P.

, and

Wood

,

R. J.

,

2010

, “

Hyperelastic Pressure Sensing With a Liquid-Embedded Elastomer

,”

J. Micromech. Microeng.

,

20

(

12

), p.

125029

.

Google Scholar

Crossref

19.

Park

,

Y.-L.

,

Tepayotl-Ramirez

,

D.

,

Wood

,

R. J.

, and

Majidi

,

C.

,

2012

, “

Influence of Cross-Sectional Geometry on the Sensitivity and Hysteresis of Liquid-Phase Electronic Pressure Sensors

,”

Appl. Phys. Lett.

,

101

(

19

), p.

191904

.

Google Scholar

Crossref

20.

Tepáyotl-Ramírez

,

D.

,

Lu

,

T.

,

Park

,

Y.-L.

, and

Majidi

,

C.

,

2013

, “

Collapse of Triangular Channels in a Soft Elastomer

,”

Appl. Phys. Lett.

,

102

(

4

), p.

044102

.

Google Scholar

Crossref

21.

Huh

,

D.

,

Mills

,

K.

,

Zhu

,

X.

,

Burns

,

M. A.

,

Thouless

,

M.

, and

Takayama

,

S.

,

2007

, “

Tuneable Elastomeric Nanochannels for Nanofluidic Manipulation

,”

Nat. Mater.

,

6

(

6

), pp.

424

–

428

.

Google Scholar

Crossref

PubMed

22.

Sparreboom

,

W.

,

Van Den Berg

,

A.

, and

Eijkel

,

J.

,

2009

, “

Principles and Applications of Nanofluidic Transport

,”

Nat. Nanotechnol.

,

4

(

11

), pp.

713

–

720

.

Google Scholar

Crossref

PubMed

23.

Ozutemiz

,

K. B.

,

Wissman

,

J.

,

Ozdoganlar

,

O. B.

, and

Majidi

,

C.

,

2018

, “

Egain–Metal Interfacing for Liquid Metal Circuitry and Microelectronics Integration

,”

Adv. Mater. Interfaces

,

5

(

10

), p.

1701596

.

Google Scholar

Crossref

24.

Sadd

,

M. H.

,

2009

,

Elasticity: Theory, Applications, and Numerics

, 2nd ed.,

Academic Press

, Waltham, MA.

25.

Muskhelishvili

,

N. I.

,

1977

,

Some Basic Problems of the Mathematical Theory of Elasticity

, 1st ed.,

Springer

, Berlin.

26.

Smirnov

,

V. I.

,

1964

,

A Course of Higher Mathematics

, Vol.

5

,

Elsevier

, Oxford, England.

27.

Savin

,

G. N.

,

1961

,

Stress Concentration Around Holes

, Vol.

1

,

Pergamon Press

, Oxford, England.

28.

Kachanov

,

M. L.

,

Shafiro

,

B.

, and

Tsukrov

,

I.

,

2003

,

Handbook of Elasticity Solutions

,

Springer Science & Business Media

, Berlin, Germany.

29.

Anderson

,

T. L.

, and

Anderson

,

T.

,

2005

,

Fracture Mechanics: Fundamentals and Applications

,

CRC Press

, Bacon Raton, FL.

30.

Adams

,

G. G.

,

2015

, “

Critical Value of the Generalized Stress Intensity Factor for a Crack Perpendicular to an Interface

,”

Proc. R. Soc. A

,

471

, p.

20150571

.

Google Scholar

Crossref

Deformation of Microchannels Embedded in an Elastic Medium

Introduction

Theory

Numerical Solution.

Analytic Solution.

Approximate Solution.

Elliptical Cross Section.

Rectangular Cross Section.

Results and Discussion

Conclusion

Acknowledgment

Funding Data

References

Get Email Alerts

Cited By

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

Deformation of Microchannels Embedded in an Elastic Medium

Introduction

Theory

Numerical Solution.

Analytic Solution.

Approximate Solution.

Elliptical Cross Section.

Rectangular Cross Section.

Results and Discussion

Conclusion

Acknowledgment

Funding Data

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only