Theory of the flow-induced deformation of shallow compliant microchannels with thick walls

TL;DR

This paper presents a theoretical model for the steady-state deformation of shallow, thick-walled microchannels under fluid flow, providing a new flow rate-pressure drop relation that aligns with experimental data.

Contribution

It develops a parameter-free analytical model for thick-walled microchannel deformation using lubrication and elasticity theories, extending previous models.

Findings

Deformation can be considered independently at each cross-section.

The top wall behaves as a simply supported rectangle under uniform pressure.

Derived flow rate-pressure drop relation matches experimental results.

Abstract

Long, shallow microchannels embedded in thick soft materials are widely used in microfluidic devices for lab-on-a-chip applications. However, the bulging effect caused by fluid--structure interactions between the internal viscous flow and the soft walls has not been completely understood. Previous models either contain a fitting parameter or are specialized to channels with plate-like walls. This work is a theoretical study of the steady-state response of a compliant microchannel with a thick wall. Using lubrication theory for low-Reynolds-number flows and the theory for linearly elastic isotropic solids, we obtain perturbative solutions for the flow and deformation. Specifically, only the channel's top wall deformation is considered, and the ratio between its thickness $t$ and width $w$ is assumed to be $(t/w)^2 \gg 1$. We show that the deformation at each stream-wise cross-section can be considered independently, and that the top wall can be regarded as a simply supported rectangle subject to uniform pressure at its bottom. The stress and displacement fields are found using Fourier series, based on which the channel shape and the hydrodynamic resistance are calculated, yielding a new flow rate--pressure drop relation without fitting parameters. Our results agree favorably with, and thus rationalize, previous experiments.

Tables

Table 1. Table 1 : Values of physical parameters used in the experiments of Gervais et al.   [ 21 ] , where the Poisson ratio is ν = 0.5 𝜈 0.5 \nu=0.5 , and the fluid’s viscosity is μ = 0.001 𝜇 0.001 \mu=0.001 Pa ⋅ ⋅ \cdot s for all the cases.

Case	$h_0$	$w$	$l$	$t$	$E_Y$	$\delta$	$ϵ$	$\gamma$	$1/{\left(\gamma \delta \right)}^2$
	[$\mu$m]	[$\mu$m]	[cm]	[mm]	[MPa]	[ – ]	[ – ]	[ – ]	[ – ]
GEGJ 1 ($▲$)	26	250	1	2	2.2	0.1040	0.0026	76.92	0.0156
GEGJ 2 ($■$)	30	500	1	2	2.2	0.0600	0.0030	66.67	0.0625
GEGJ 3 ($▲$)	26	250	1	2	1.1	0.1040	0.0026	76.92	0.0156
GEGJ 4 ($■$)	30	500	1	2	1.1	0.0600	0.0030	66.67	0.0625

Table 2. Table 2 : The comparison of the values of R 2 superscript 𝑅 2 R^{2} between the present theory, equation ( 28 ), and the model in [ 21 ] , denoted GEGJ fit (see equation ( 29 )). The value of the fitting parameter α 𝛼 \alpha for each case is available in [ 21 ] .

Case	GEGJ 1	GEGJ 2	GEGJ 3	GEGJ 4
Present theory	0.9633	0.9859	0.8904	0.9224
GEGJ fit	0.9904	0.9988	0.9792	0.9920

Table 3. Table 3 : Values of the physical parameters used in the experiments of Raj et al.   [ 26 ] . The microchannel is w = 350 𝑤 350 w=350 μ 𝜇 \mu m wide, l = 3 𝑙 3 l=3 cm long and h 0 = 50 subscript ℎ 0 50 h_{0}=50 ± plus-or-minus \pm 5 μ 𝜇 \mu m in height. Based on the reported experimental conditions in [ 26 ] , the fluid viscosity is taken to be μ = 9.110 × 10 − 4 𝜇 9.110 superscript 10 4 \mu=9.110\times 10^{-4} Pa ⋅ ⋅ \cdot s, and the Poisson ratio is ν = 0.5 𝜈 0.5 \nu=0.5 , for all the cases.

Case	$t$	$E_Y$	$\delta$	$ϵ$	$\gamma$	$1/{\left(\gamma \delta \right)}^2$
	[mm]	[MPa]	[ – ]	[ – ]	[ – ]	[ – ]
RDC 1 ($\square$)	2.0	2.801	0.1286–0.1571	0.0015–0.0018	36.36–44.44	0.0306
RDC 2 ($\square$)	1.0	2.801	0.1286–0.1571	0.0015–0.0018	18.18–22.22	0.1225
RDC 3 ($\square$)	0.5	2.801	0.1286–0.1571	0.0015–0.0018	9.091–11.11	0.49
RDC 4 ($○$)	2.0	0.157	0.1286–0.1571	0.0015–0.0018	36.36–44.44	0.0306
RDC 5 ($○$)	1.0	0.157	0.1286–0.1571	0.0015–0.0018	18.18–22.22	0.1225
RDC 6 ($○$)	0.5	0.157	0.1286–0.1571	0.0015–0.0018	9.091–11.11	0.49

Equations

$$X=\frac{x}{w},Y=\frac{y}{h_0},Z=\frac{z}{l},V_X=\frac{\delta v_x}{ϵ{\mathcal{V}}_c},V_Y=\frac{v_y}{ϵ{\mathcal{V}}_c},V_Z=\frac{v_z}{{\mathcal{V}}_c},P=\frac{p}{{\mathcal{P}}_0},$$

$$\frac{\partial V_X}{\partial X}+\frac{\partial V_Y}{\partial Y}+\frac{\partial V_Z}{\partial Z}$$

$$\frac{{ϵ}^3}{{\delta }^2}Re\left(V_X\frac{\partial V_X}{\partial X}+V_Y\frac{\partial V_X}{\partial Y}+V_Z\frac{\partial V_X}{\partial Z}\right)$$

$${ϵ}^{3}Re\left(V_X\frac{\partial V_Y}{\partial X}+V_Y\frac{\partial V_Y}{\partial Y}+V_Z\frac{\partial V_Y}{\partial Z}\right)$$

$$ϵRe\left(V_X\frac{\partial V_Z}{\partial X}+V_Y\frac{\partial V_Z}{\partial Y}+V_Z\frac{\partial V_Z}{\partial Z}\right)$$

$$V_Z(X,Y,Z)=\frac{1}{2}\frac{\mathrm{d}P}{\mathrm{d}Z}(Y+1)[Y-H(X,Z)](-1\le Y\le H).$$

$$v_z(x,y,z)=\frac{1}{2\mu }\frac{\mathrm{d}p}{\mathrm{d}z}(y+h_0)\left[y-u_y^0(x,z)\right]\left(-h_0\le y\le u_y^0\right).$$

$$-\frac{{\mathcal{P}}_0}{l}\frac{\partial P}{\partial Z}+\frac{{\mathfrak{T}}_0}{h_0}\frac{\partial {\mathcal{T}}_{YZ}}{\partial Y}=0.$$

$$\frac{\partial {\sigma }_{xx}}{\partial x}+\frac{\partial {\sigma }_{xy}}{\partial y}+\frac{\partial {\sigma }_{xz}}{\partial z}$$

$$\frac{\partial {\sigma }_{xy}}{\partial x}+\frac{\partial {\sigma }_{yy}}{\partial y}+\frac{\partial {\sigma }_{yz}}{\partial z}$$

$$\frac{\partial {\sigma }_{xz}}{\partial x}+\frac{\partial {\sigma }_{yz}}{\partial y}+\frac{\partial {\sigma }_{zz}}{\partial z}$$

$${\Sigma }_{XX}=\frac{{\sigma }_{xx}}{{\mathcal{D}}_{xx}},{\Sigma }_{\mathcal{YY}}=\frac{{\sigma }_{yy}}{{\mathcal{P}}_0},{\Sigma }_{ZZ}=\frac{{\sigma }_{zz}}{{\mathcal{D}}_{zz}},{\Sigma }_{X\mathcal{Y}}=\frac{{\sigma }_{xy}}{{\mathcal{D}}_{xy}},{\Sigma }_{\mathcal{Y}Z}=\frac{{\sigma }_{yz}}{ϵ{\mathcal{P}}_0}.$$

$$\frac{{\mathcal{D}}_{xx}}{w}\frac{\partial {\Sigma }_{XX}}{\partial X}+\frac{{\mathcal{D}}_{xy}}{\gamma h_0}\frac{\partial {\Sigma }_{X\mathcal{Y}}}{\partial \mathcal{Y}}$$

$$\frac{{\mathcal{D}}_{xy}}{w}\frac{\partial {\Sigma }_{X\mathcal{Y}}}{\partial X}+\frac{{\mathcal{P}}_0}{\gamma h_0}\frac{\partial {\Sigma }_{\mathcal{YY}}}{\partial \mathcal{Y}}+\frac{ϵ{\mathcal{P}}_0}{l}\frac{\partial {\Sigma }_{\mathcal{Y}Z}}{\partial Z}$$

$$\frac{ϵ{\mathcal{P}}_0}{\gamma h_0}\frac{\partial {\Sigma }_{\mathcal{Y}Z}}{\partial \mathcal{Y}}+\frac{{\mathcal{D}}_{zz}}{l}\frac{\partial {\Sigma }_{ZZ}}{\partial Z}$$

$$\frac{{\partial }^4\varphi }{\partial x^4}+2\frac{{\partial }^4\varphi }{\partial x^2\partial y^2}+\frac{{\partial }^4\varphi }{\partial y^4}=0.$$

$${\sigma }_{yy}{|}_{y=0}=-p(z)=-p(z)\sum _{m=1}^{\infty }{A}_m\sin \left[{\kappa }_m\left(x+\frac{w}{2}\right)\right],$$

$${\varphi }_m(x,y)=\sin \left[{\kappa }_m\left(x+\frac{w}{2}\right)\right]\left(C_1{e}^{{\kappa }_{m}y}+C_2{e}^{-{\kappa }_{m}y}+C_{3}y{e}^{{\kappa }_{m}y}+C_{4}y{e}^{-{\kappa }_{m}y}\right)$$

$${\sigma }_{xx,m}(x,y)$$

$$\displaystyle\quad+C_{3}(2\kappa_{m}e^{\kappa_{m}y}+e^{\kappa_{m}y})+C_{4}(-2\kappa_{m}e^{-\kappa_{m}y}+\kappa_{m}^{2}ye^{-\kappa_{m}y})\big{]},$$

$${\sigma }_{yy,m}(x,y)$$

$${\sigma }_{xy,m}(x,y)$$

$$\displaystyle\quad+C_{3}(e^{\kappa_{m}y}+{\kappa_{m}}ye^{\kappa_{m}y})+C_{4}(e^{-\kappa_{m}y}-{\kappa_{m}}ye^{-\kappa_{m}y})\big{]}.$$

$${\sigma }_{yy,m}{|}_{y=0}=A_m\sin \left[{\kappa }_m\left(x+\frac{w}{2}\right)\right],{\sigma }_{xy,m}{|}_{y=0}=0,{\sigma }_{yy,m}{|}_{y=t}=0,{\sigma }_{xy,m}{|}_{y=t}=0.$$

$$C_1$$

$$C_2$$

$$C_3$$

$$C_4$$

$$\begin{gathered} \varphi (x,y,z)=-p(z)\sum _{m=1}^{\infty }{\varphi }_m(x,y),{\sigma }_{xx}(x,y,z)=-p(z)\sum _{m=1}^{\infty }{\sigma }_{xx,m}(x,y), \\ {\sigma }_{yy}(x,y,z)=-p(z)\sum _{m=1}^{\infty }{\sigma }_{yy,m}(x,y),{\sigma }_{xy}(x,y,z)=-p(z)\sum _{m=1}^{\infty }{\sigma }_{xy,m}(x,y). \end{gathered}$$

$$e_{xx,m}$$

$$e_{yy,m}$$

$$e_{xy,m}$$

$$u_{x,m}(x,y)$$

$$\displaystyle\phantom{=}+[2C_{3}+\kappa_{m}(1+\overline{\nu})(C_{1}+C_{3}y)]e^{\kappa_{m}y}\Big{\}}+f_{2}(y),$$

$$u_{y,m}(x,y)$$

$$\displaystyle\phantom{=}-[-C_{3}(1-\overline{\nu})+\kappa_{m}(1+\overline{\nu})(C_{1}+C_{3}y)]e^{\kappa_{m}y}\Big{\}}+f_{1}(x),$$

$$\frac{\partial f_1}{\partial x}+\frac{\partial f_2}{\partial y}=0.$$

$$u_x(x,y,z)=-p(z)\sum _{m=1}^{\infty }{u}_{x,m}(x,y),u_y(x,y,z)=-p(z)\sum _{m=1}^{\infty }{u}_{y,m}(x,y).$$

$$u_y^0(x,z)=\frac{p(z)}{{\overline{E}}_Y}\sum _{m=1}^{\infty }\frac{2A_m}{{\kappa }_m}\sin \left[{\kappa }_m\left(x+\frac{w}{2}\right)\right]({\gamma }^2{\delta }^2=t^2/w^2\to \infty ).$$

$$\frac{u_c{\overline{E}}_Y}{w{\mathcal{P}}_0}{U}_Y^0(X,Z)=P(Z)\sum _{m=1}^{\infty }\frac{2A_m}{m\pi }\sin \left[m\pi \left(X+\frac{1}{2}\right)\right]=P(Z)\mathfrak{G}(X),$$

Full text

{fmtext}

1 Introduction

Microfluidic devices have enabled the miniaturization of processes that involve the flow and manipulation of small volumes of fluids, down to the nanoliter [1]. In recent years, polymer-based materials have become popular for the fabrication of microfluidic devices because they promise cheaper and faster production cycles [2]. A widely-used material for microfabrication is

the elastomer poly(dimethylsiloxane) (PDMS) [3, 4, 5]. The development of PDMS-based microfluidic devices is intimately related to the emergence of the cutting-edge field known as lab-on-a-chip [6, 7]. Microchannels made from PDMS have found many applications, such as platforms for organ-on-a-chip models [8] and various biological studies (e.g., assays and combinatorial screening) [9]. PDMS, also known commercially as Sylgard 184, typically has a low tensile modulus [10], and PDMS-based micro-conduits are prone to deformation (“bulging”) under applied forces. Although early studies considered deformation to be a drawback, because it may restrict the structural viability of a device [11, 12], the compliance of PDMS microchannels has been exploited to design microfluidic devices with specific functions, such as pressure-actuated valves [13], passive fuses [14], pressure sensors [15, 16], strain sensors [17], and impedance-based flow meters with improved sensitivity [18].

The soft nature of microfluidic devices made from PDMS or similar materials requires the consideration of fluid–structure interactions (FSIs), even at the low Reynolds numbers encountered at these scales [19]. Given the broad range of applications of these microsystems, understanding the fluid–structure interplay is important and necessary [20]. Conceptually speaking, the soft wall of the microchannel will bulge due to the hydrodynamic pressure within, and this effect will, in turn, modify the flow velocity and pressure gradient. Rectangular microchannels embedded in soft materials have been the most common configuration. However, even within this “simple” Cartesian geometry, quantifying the FSI loop is not easy. To this end, in the present work, we analyze the steady-state FSI in a long and shallow compliant microchannel with thick walls. Gervais et al. [21] were the first to put forward a model for this case by introducing a fitting parameter to be determined via calibration with experiments. In their model, the strain in the elastic solid is taken to be linearly proportional to the hydrodynamic pressure, with the proportionality constant being unknown. The same scaling analysis was generalized to the unsteady case by Dendukuri et al. [22] for studying stop-flow lithography in a thick-walled microfluidic device. Mukherjee et al. [23] then considered actuation of the soft wall via electroosmotic flow, eliminating the fitting parameter under the assumption of two-dimensional (2D) kinematics (i.e., no displacement in the spanwise direction transverse to the flow). Without the assumption of 2D kinematics, the cross-sectional deformation profile in the spanwise direction was previously assumed to be parabolic. In this way, these studies were able to quantify the flow rate–pressure drop relation and show that this relation deviates from the linear proportionality predicted by Poiseuille’s law. Although the model in [21] was initially developed for microchannels with walls behaving like a semi-infinite medium, the model has actually (with varying degrees of success) been also applied to FSI in microchannels with elastic walls of various thicknesses [24, 25, 26].

The fitting parameter in the model from [21] is inconvenient in applications because it has to be recalibrated for each microchannel’s geometric and material properties. This model also neglects many fluid–solid coupling details because a solution of the elasticity problem has not been obtained from the governing equations. It has been previously noticed that the fitting parameter has to be related to the top wall thickness [24]. Seker et al. [27] argued that half-space-like thick walls and plate-like thin walls should be treated differently, and they used a empirical formula to determine the fitting parameter in these two regimes. Later, Christov et al. [28] used perturbation methods and successfully obtained a fitting-parameter-free model for the flow rate–pressure drop relation in a microchannel with a clamped plate-like top wall. This work was then extended by Shidhore and Christov [29] to account for the deformation of microchannels whose top wall behaves like a thick plate. Most recently, Martínez-Calvo et al. [30] considered the unsteady version of the problem in [28] in order to update the previous fitting-parameter-based approach [22]. This line of research has shown that fitting parameters are not necessary to quantify microscale FSIs, and they can be avoided entirely by seeking a perturbative solution to the related elasticity problem. All of these models mentioned are applicable to Newtonian fluid flow only. Other works have analyzed FSIs in microchannels conveying non-Newtonian fluids [31, 32].

All of the existing fitting-parameter-free models target microchannels with thin walls (i.e., plate-like elastic structures) and, thus, fail immediately when pushed to a large-thickness regime. Yet, standard fabrication techniques based on soft lithography [3, 5] naturally produce microchannels with thick elastic slabs as the boundary, violating both the plate-like assumption and the clamped-boundary assumption. The present work seeks to fill this knowledge gap between theory and experiment. To this end, we derive a new mathematical model (without fitting parameters) for the steady fluid–structure interaction between a Newtonian fluid and an (initially) rectangular microchannel with deformable walls. Specifically, via a scaling analysis of the equilibrium stress balance in the elastic solid, we give a precise definition of “thick-walled” microchannels in an asymptotic sense (section 3.1). The reduced solid mechanics problem is then solved exactly by Fourier series (section 3.2), yielding the precise shape of the fluid–solid interface profile (section 3.3). Next, the flow rate–pressure drop relation is calculated (section 4) by substituting the deformation profile into the leading-order flow field under the lubrication approximation (section 2). The final mathematical result for the nonlinear hydrodynamic resistance is compared with experimental measurements reported in [21, 26], and it is found to be in favorable agreement (section 5). Finally, conclusions and avenues for future work are discussed in section 6. For completeness, Appendix is provided discussing how the thick-slab elasticity problem limits onto the plate-like problem from previous theoretical studies.

Problem statement, geometry and notation

Consider the pressure-driven, steady flow of a Newtonian viscous fluid inside of a microchannel fabricated via soft lithographic techniques from an elastomer. A schematic configuration is shown in figure 1. The problem is symmetric about $x=0$, with only half of it being shown in the figure. The fluid domain within microchannel has undeformed height $h_{0}$, width $w$ and length $l$. The flow is in the positive $z$-direction. It is assumed that the microchannel is long and shallow, such that $h_{0}\ll w\ll l$. Three sides (two lateral walls and the top wall) of the channel are composed of a soft elastic material, while the bottom wall is assumed to be rigid. Here, $t$ denotes the thickness of the top solid slab, while the thickness of the lateral walls is assumed to be large enough to be considered as infinite. According to the scaling analysis of Gervais et al. [21], the strain of the displaced fluid–solid interface is proportional to the imposed stress, so $u_{y}^{0}(x,z)\propto wp(z)/E_{Y}$, where $p(z)$ denotes the local pressure at flow-wise position $z$, and $E_{Y}$ denotes the Young’s modulus of the solid. Similarly, the side wall deformation is proportional to $h_{0}p(z)/E_{Y}$. In the regime $h_{0}\ll w$, the deformation of the top wall is thus expected to be much larger than that of the lateral walls. Therefore, we only consider the deformation of the top wall.

2 Governing equations of the fluid mechanics problem

First, let us introduce the dimensionless variables

[TABLE]

where $\epsilon\mathcal{V}_{c}/\delta$, $\epsilon\mathcal{V}_{c}$ and $\mathcal{V}_{c}$ are respectively the characteristic scales for the velocity components $v_{x}$, $v_{y}$ and $v_{z}$, chosen to balance the conservation of mass equation (2a). Meanwhile, $\mathcal{P}_{0}$ is the characteristic pressure scale. If the pressure is controlled at both the inlet and outlet, then $\mathcal{P}_{0}=\Delta p=p(z=l)-p(z=0)$, which is the total pressure drop over the channel length $l$, and $\mathcal{V}_{c}=h_{0}^{2}\Delta p/(\mu l)$. Instead, if the volumetric flow rate $q$ is controlled at the inlet (and the flow is steady, as assumed), then $\mathcal{V}_{c}$ is chosen as the mean flow velocity, $\langle v_{z}\rangle=q/(h_{0}w)$, and $\mathcal{P}_{0}=\mu\langle v_{z}\rangle l/h_{0}^{2}$; see the detailed discussion in [28]. We have also defined the dimensionless ratios $\epsilon=h_{0}/l$, $\delta={h_{0}}/w$, with $\epsilon\ll\delta\ll 1$ (recall section 1), and the Reynolds number $Re=\rho\mathcal{V}_{c}h_{0}/{\mu}$, where $\rho$ and $\mu$ are the fluid density and dynamic viscosity, respectively. Assuming the working fluid is incompressible, Newtonian, and steady, then the governing equations are the steady incompressible Navier–Stokes equations, which take the dimensionless form:

[TABLE]

As shown in figure 1, the fluid mechanics problem is posed on the deformed domain $\Omega_{f}=\{(x,y,z)\,|-w/2\leq x\leq+w/2,\,-h_{0}\leq y\leq u_{y}^{0}(x,z),\,0\leq z\leq l\}$. We set $u_{y}^{0}=u_{c}U_{Y}^{0}$ and $\lambda=u_{c}/h_{0}$, where $u_{c}$ is the characteristic top wall deformation scale (to be determined self-consistently by solving the corresponding elasticity problem in section 3). Then, $H(X,Z)=\lambda U_{Y}^{0}(X,Z)$ is the dimensionless deformed top wall profile, and $\lambda$ can be interpreted as the compliance parameter that characterizes top wall’s ability to deform due to the flow beneath it. Under the lubrication approximation, which we shall now introduce, it is expected that $\lambda\ll 1/\delta$ [28]. Further assuming that $\epsilon Re\ll 1$, the leading-order solution of equations (2), subject to no-slip condition at the top and bottom walls of the channel, is

[TABLE]

At this stage, the (non-constant) pressure gradient $\mathrm{d}P/\mathrm{d}Z<0$ remains unknown. Due to the chosen scaling (to balance conservation of mass), the velocity components $V_{X}$ and $V_{Y}$ come in at higher orders in the perturbation expansion. Systematic corrections in powers of $\epsilon$ can be obtained as regular perturbations [33], however the expansion in $\delta$ is singular [34]. For our purposes, it suffices to note that the flow is primarily unidirectional in the $z$-direction. This is the familiar Reynolds lubrication approximation [35, 36]. It is also useful to write down the dimensional form of equation (3) for the convenience in the upcoming discussion:

[TABLE]

At the leading order, the balance between the fluid’s shear stress and the pressure ensures the conservation of linear momentum in the $z$-direction (see equation (2d)). It is useful to rewrite the equation explicitly in the form of stresses and make dimensionless in order to determine the ordering of the forces, which reads

[TABLE]

Here, we have introduced $\mathfrak{T}_{0}$ as the characteristic scale for the fluid’s shear stress $\mathcal{T}_{YZ}$. The force balancing requires the prefactors to be of the same order, thus $\mathfrak{T}_{0}=\epsilon\mathcal{P}_{0}$, i.e., the shear stress in the flow is much smaller than the pressure. Therefore, when dealing with the solid mechanics problem in the next section, the hydrodynamic pressure is regarded as the only applied force that causes deformation.

3 Governing equations of the solid mechanics problem

3.1 Plane strain configuration and the thickness effect

It should be clarified beforehand that we “cut off” the solid from the sides and only consider the deformation of the top rectangular slab, with width $w$, thickness $t$ and length $l$ with initial configuration occupying the domain $\Omega_{s0}=\{(x,y,z)\,|-w/2\leq x\leq+w/2,\,0\leq y\leq t,\,0\leq z\leq l\}$. For convenience, we introduce the parameter $\gamma=t/h_{0}$ to denote the ratio of the solid thickness to the undeformed channel height and assume that $t\ll l$, equivalent to $\gamma\epsilon\ll 1$. Furthermore, the deformation is assumed to be much smaller than any dimension of the solid (i.e., small deformation gradient) so that the theory of linear elasticity applies, i.e., $u_{c}\ll w$ and $u_{c}\ll t$. In terms of the compliance parameter $\lambda$, this assumption requires that $\lambda\ll\mathrm{max}(1/\delta,\gamma)$. Neglecting body forces, the elastostatics equations are simply reduced to the balance between the Cauchy stresses:

[TABLE]

where $\sigma_{xx}$, $\sigma_{yy}$, $\sigma_{zz}$, $\sigma_{xy}$, $\sigma_{xz}$ and $\sigma_{yz}$ are the six independent variables.

The next natural step is to find the primary force balancing by making the elastostatics equations dimensionless. First, the scale for $y$ should now be taken as $\gamma h_{0}=t$ in order to account for the thickness effect, thus let $\mathcal{Y}=y/t$. Second, the continuity of stresses at the fluid–solid interface gives a clue of how to pick the scales for the stresses in the solid. It has been shown in section 2 that $p\sim\mathcal{P}_{0}$ and $\tau_{yz}\sim\epsilon\mathcal{P}_{0}$. Since $\sigma_{yy}=-p$ and $\tau_{yz}=\sigma_{yz}$ at the fluid–solid interface, the scales for $\sigma_{yy}$ and $\sigma_{yz}$ should be $\mathcal{P}_{0}$ and $\epsilon\mathcal{P}_{0}$, respectively. What is more, it is believed that $\sigma_{xz}$ is negligible because $\tau_{xz}=\mathcal{O}(\epsilon\delta)$ is even smaller. Thus, we henceforth neglect $\sigma_{xz}$ from the equilibrium equations (6). Accordingly, we make the stress components of the solid dimensionless as

[TABLE]

where $\mathcal{D}_{xx}$, $\mathcal{D}_{zz}$, and $\mathcal{D}_{xy}$ are the (a priori unknown) characteristic scales for $\sigma_{xx}$, $\sigma_{zz}$ and $\sigma_{xy}$, respectively. Then the dimensionless elastostatics equations become

[TABLE]

From equation (8b), the normal stress has to be balanced with the shear stress, at the leading order in $\epsilon$, thus $\mathcal{D}_{xy}=\mathcal{P}_{0}/(\gamma\delta)$. Substituting this scale into equations (8a) and (8c), we obtain $\mathcal{D}_{xx}=\mathcal{P}_{0}/(\gamma^{2}\delta^{2})$ and $\mathcal{D}_{zz}=\mathcal{P}_{0}/\gamma$. It is important to note that $\gamma$ (i.e., the solid thickness parameter) plays an essential role in the stress distribution. Accordingly, the boundary condition at the sidewalls $x=\pm w/2$, due to the reaction between the top wall and the remaining solid, will give rise to the thickness effect.

In the present study, we are interested in the thickness range of $t\ll l$, i.e., $\gamma\ll 1/\epsilon$, thus $\gamma\delta\ll\delta/\epsilon$. In this distinguished limit, we have shown that $\sigma_{xx}$, $\sigma_{xy}$, $\sigma_{yy}$ and $\sigma_{zz}$ are the dominant stresses in the solid as $\epsilon\to 0^{+}$. Furthermore, based on the linear constitutive relation between the stress and the linear strain, as well as the fact that the microchannel is usually prevented from displacements in the flow-wise direction by rigid inlet and outlet connectors, we are justified in reducing the problem to a plane strain configuration with dominant linear strains $e_{xx}$, $e_{yy}$ and $e_{yz}$, in the cross-section.

For a plane strain problem, it is convenient to introduce the Airy stress function $\phi(x,y)$ that satisfies the homogeneous biharmonic equation (in dimensional form) [37]:

[TABLE]

Then, the stresses are computed from the stress function as $\sigma_{xx}=\partial^{2}\phi/\partial y^{2}$, $\sigma_{yy}=\partial^{2}\phi/\partial x^{2}$ and $\sigma_{xy}=-\partial^{2}\phi/\partial x\partial y$. Meanwhile, $\sigma_{zz}=\nu(\sigma_{xx}+\sigma_{yy})$ according to the constitutive equation of linear elasticity, where $\nu$ is the Poisson ratio. Note our analysis necessitates different characteristic scales for the $y$ coordinate for the fluid and solid mechanics problems. Therefore, for consistency and convenience, we solve the solid mechanics problem in its dimensional form (9).

3.2 Large-thickness case

In the present study, as discussed in section 1, we are interested in the case of thick top wall, which arises because microchannels are frequently embedded in half-space-like PDMS medium when manufactured by, e.g., replica molding [5]. From the scaling analysis of equation (8) above, we learned that, as the thickness increases, $\sigma_{xy}$ decreases as $1/(\gamma\delta)$, while $\sigma_{xx}$ decreases even faster, as $1/(\gamma\delta)^{2}$. Thus, consider the case when the thickness is large enough, specifically $(\gamma\delta)^{2}\gg 1$. Then, $\sigma_{xx}$ is much smaller than the other stresses in the solid, as well as at both sidewalls. Taking $\sigma_{xx}|_{x=\pm w/2}=0$ and assuming that the displacement at the corner is negligible, the boundary condition at $x=\pm w/2$ is reduced to that of a simple support. This result is crucial to the analysis that follows because equation (9) can now be solved exactly using Fourier series in the simply supported rectangular configuration.

First, the stress at the fluid–solid interface is decomposed into Fourier series. At $y=0$, the normal stress in the solid should match the local hydrodynamic pressure:

[TABLE]

where $\kappa_{m}={m\pi}/{w}$ and $A_{m}=\frac{2}{m\pi}[1-(-1)^{m}]$. We have used the fact that the Fourier series part of (10) equals $1$ for $x\in(-w/2,+w/2)$. Note, however, that the series converges to [math] at $x=\pm{w}/{2}$ because the odd extension has been used here to construct sine series, which causes the discontinuity at the two edges.

Next, the superposition principle comes into play. It is easily verified that

[TABLE]

satisfies equation (9) for any integer $m=1,2,\ldots$. The corresponding stress state is

[TABLE]

Four boundary conditions are needed to determine these coefficients. The stress continuity at the interface (equations (13)1 and (13)2) and the stress free conditions at the upper edge of the top wall (equations (13)3 and (13)4) require that

[TABLE]

Imposing equations (13) on equations (12), we obtain

[TABLE]

Note the coefficients $C_{1}$, $C_{2}$, $C_{3}$ and $C_{4}$ are not fixed constants but vary with $m$ and the thickness $t$ via the definition $\beta=\kappa_{m}t=m\pi\gamma\delta$. Finally, the solution to equation (9), as well as the three unique stress components, can be constructed by superposition:

[TABLE]

3.3 Displacements and the shape of the fluid–solid interface

Of course, the analysis above is only valid for small deformation gradients. In this regime, the stress–strain relations of linear elasticity [37] dictate that

[TABLE]

where $G=E/[2(1+\nu)]$ is the shear modulus of elasticity, and $\overline{E}_{Y}$ and $\overline{\nu}$ are related to the Young’s modulus $E_{Y}$ and the Poisson’s ratio $\nu$ by $\overline{E}_{Y}={E_{Y}}/(1-\nu^{2})$ and $\overline{\nu}={\nu}/(1-\nu)$, respectively, because of the plane strain configuration considered herein [37].

Integrating equations (16a) and (16b), $u_{x,m}$ and $u_{y,m}$ are, respectively,

[TABLE]

where $f_{1}(x)$ and $f_{2}(y)$ are arbitrary functions of integration. Substituting equations (17) and (17) into equation (16c), we find

[TABLE]

From equation (18), it is easily concluded that both $f_{1}$ and $f_{2}$ should be constants. Since it is assumed that there are no vertical displacement at $x=\pm{w}/{2}$, $f_{1}=0$. According to the symmetry of the problem, i.e., $({\partial u_{x}}/{\partial x})|_{x=0}=0$, $f_{2}=0$ as well.

Finally, the displacements are obtained by summing up all the $u_{x,m}$ and $u_{y,m}$ terms from equations (17):

[TABLE]

To obtain the fluid–solid interface deflection profile, $u_{y}^{0}(x,y)$, as well as the shape of the whole deformed cross-section, we simply take $y=0$ in equations (19). In order to evaluate the Fourier series numerically and generate the plots herein, we find that keeping 50 terms in the sum is sufficient.

However, in the large thickness case of interest herein, the above results can be further simplified because $C_{1}$ and $C_{3}$ are small compared to $C_{2}$ and $C_{4}$, respectively. As shown in figure 2, the dimensionless ratios $C_{1}/C_{2}$ and $C_{3}/C_{4}$ decrease very quickly with $\gamma\delta$. Specifically, for $t/w\gtrsim 1$ (say, $t\simeq 1.5w$), $C_{1}$ and $C_{3}$ are negligible compared to $C_{2}$ and $C_{4}$. In this case, we can simply regard the stress-free boundary conditions in equations (13)3 and (13)4 as being satisfied at $y=\infty$ instead of $y=t$; hence, $C_{1}=C_{3}=0$, $C_{2}=-A_{m}/\kappa_{m}^{2}$ and $C_{4}=-A_{m}/\kappa_{m}$. Then, the vertical displacement at the fluid–solid interface reduces to

[TABLE]

The panel of figure 2 supports our observations. The interface profiles for the cases of $\gamma\delta=1.5$ and $\gamma\delta=2.0$ coincide with the curve predicted by equation (20).

It is easy to rewrite equation (20) in dimensionless form as

[TABLE]

where, for convenience, we have denoted by $\mathfrak{G}$ the function of $X$ defined by the Fourier series. Now, the natural deformation scale is clearly $u_{c}=w\mathcal{P}_{0}/\overline{E}_{Y}$, so that we can set the prefactor on the left-hand side of equation (21) to unity. This scale is similar to the one used in [21], where it was assumed $\langle u_{y}\rangle/w\sim p/E_{Y}$. Note, however, that our analysis shows that $\overline{E}_{Y}=E_{Y}/(1-\nu^{2})$ must be used in the deformation scale instead of $E_{Y}$ because the top wall in a long, shallow microchannel is in a plane strain configuration. Then, $\lambda$ in equation (3) is finally determined to be $\lambda=u_{c}/h_{0}=\mathcal{P}_{0}/(\overline{E}_{Y}\delta)$.

We can see from equation (21) that the interface deflection profiles at different $Z$ coordinates have the same shape, denoted by $\mathfrak{G}(X)$. It is easy to compute the maximum and average displacement at the interface from equation (21):

[TABLE]

Note the prefactors here are different from previous studies, which either assumed a parabolic deformation profile of the fluid–solid interface, in which case $\langle U_{Y}^{0}\rangle=(2/3)\max_{X}U_{Y}^{0}$ [16], or obtained a quartic profile from plate theory with clamped boundary condition, in which case $\langle U_{Y}^{0}\rangle=(8/15)\max_{X}U_{Y}^{0}$ [28].

Observe that the simple support does not restrict the horizontal displacement (see equation (19)1). Denoting the horizontal displacement at the fluid–solid interface by $u_{x}^{0}$, we can express this displacement in the large-thickness case as

[TABLE]

Given that the typical Poisson ratio of PDMS is $\nu\approx 0.4$ to 0.5 [10], our theory predicts $u_{x}^{0}\ll u_{y}^{0}$ since the ratio of the maximum value of these two displacements, from equations (21) and (23), is only $(1-\overline{\nu})/2=(1-2\nu)/[2(1-\nu)]$. This result is consistent with experimental observations [21]. Interestingly, if we take the material of the top wall to be strictly incompressible with $\nu=0.5$, the $u_{x}^{0}$ is exactly zero, even though $u_{x}\neq 0$ for $y>0$. The latter is not important in the context of the present study because we focus on the fluid domain’s shape inside the microchannel.

3.4 Summary and discussion of the solid mechanics results

To summarize, we have derived a mathematical expression for the fluid–solid interface deflection curve for the large thickness case. It should be clarified again that the thickness is considered “large” specifically when $1/(\gamma\delta)^{2}=(w/t)^{2}\ll 1$. In this distinguished limit, we have shown that the top wall can be considered as a simply supported rectangle subject to uniform pressure at the bottom. Importantly, note the present thickness range includes but is wider than $w/t\ll 1$. For some cases with $w/t\simeq 1$, asymptotically, we can still satisfy $(w/t)^{2}\ll 1$.

An important case is the static response of microchannels with thinner top wall, such as a plate-like structure. However, our configuration, based on the fabrication methods from, e.g., [21] and depicted in figure 1, is different from previous studies [39, 16, 28, 29, 40, 41] in which the top wall is modelled as a clamped membrane or plate with small thickness $t\lesssim w$. Consequently, $1/(\gamma\delta)^{2}\ll 1$ is not satisfied and the analysis above is inapplicable. Nevertheless, the plane strain assumption is still valid. After analyzing the resultant forces and moments in the solid, we hypothesize that in this case, at each fixed-$z$ cross-section, the top wall behaves like a simply supported beam with tension. However, it is hard to find the exact solutions for the biharmonic equation (9) under such boundary conditions. Therefore, we use beam theory to obtain the fluid–solid interface profile and the hydrodynamic resistance. These details are provided in the Appendix.

4 Hydrodynamic resistance of the compliant channel

Having solved for the leading-order velocity profile in section 2 and the cross-sectional shape of the fluid–solid interface in section 3.3, we are now in a position to solve the coupled fluid–structure interaction problem. Specifically, in the microfluidics context, of greatest interest is the hydrodynamic resistance, which characterizes the required pressure drop (i.e., force) to maintain steady flow at a given volumetric flow rate [2]. For fixed cross-sectional shapes, this quantity can be characterized for any number of shapes [42] using the ability to solve the Stokes equations for $Re=0$ [43, section 2-5].

In a compliant channel, however, and the specific shape of the deformation profile in cross-section depends on the pressure itself. This results in a nonlinear relationship between the pressure drop and flow rate. We can determine this relation by directly calculating the flow rate under the deformed cross-section. On using equation (4) for $v_{z}(x,y,z)$, we obtain:

[TABLE]

where $u_{y}^{0}(x,z)$ is given in equation (20) for the large-thickness case, and only the axial velocity component contributes to the flow rate at the leading order in the assumed small parameters. For a scenario with constant flow rate, equation (24) is a first-order differential equation for $p(z)$ given $q=const.$, which can be solved by assuming the outlet pressure sets the gauge, i.e., $p(l)=0$.

In the large-thickness case, the self-similarity of the fluid–solid interface deflection profile makes it easy to solve equation (24) by separation of variables, yielding an implicit relation for the hydrodynamic pressure:

[TABLE]

where

[TABLE]

and $\mathfrak{G}(X)$ is the self-similar deflection profile shared by every cross-section given in equation (21). The integrals in equations (26) are computed numerically using the trapezoidal with respect to 100 evenly-space integration points on $X\in[-1/2,+1/2]$. Observe that, while equation (25) has the same general structure (as already expected from [44]) as that arising from theories based on plate-like elastic top walls [28], the prefactors $S_{1,2,3}$ related to the cross-sectional shape of the fluid–solid interface are larger by orders of magnitude.

Equation (25) can also be made dimensionless in the flow-rate-controlled regime as

[TABLE]

Recall that we defined $\lambda=u_{c}/h_{0}=\mathcal{P}_{0}/(\overline{E}_{Y}\delta)$ so that $H(X,Z)=\lambda U_{Y}^{0}(X,Z)$. The dimensionless flow rate–total pressure drop relation is obtained by taking $Z=0$ in equation (27). As the top wall deformation increases, i.e., for larger $\lambda$, the nonlinearity in the relation becomes more pronounced.

Finally, taking $z=0$, the relation between the total pressure drop $\Delta p$ and the volumetric flow rate $q$ is obtained from (25):

[TABLE]

The important message is that with the consideration of the fluid–structure interaction in the microchannel, the flow rate and pressure drop relation deviates from the classic Poiseuille’s law, which for a rectangular channel is $q=wh_{0}^{3}\Delta p/(12\mu l)$ (neglecting drag from the lateral sidewalls) [2], and displays nonlinearity. Although $S_{1}>S_{2}>S_{3}$, it is important to emphasize that equation (28) is simply a polynomial, and not a perturbation series, in $\Delta p$.

5 Results, discussion and illustrated examples

At the leading order in the small parameter $\epsilon$, we have reduced the 3D deformation of the top wall to a 2D problem by considering each cross-section (in the $(x,y)$ plane) as independent. Such decoupling is a natural consequence of the long and slender geometry and has also been shown asymptotically by Christov et al. [28] and numerically by Chakraborty et al. [45]. Based on this idea, various models, either with or without fitting parameters, have been put forward to account for the nonlinear flow rate–pressure drop relation. The very first one was from Gervais et al. [21] in the form of an implicit relation for the hydrodynamic pressure $p(z)$:

[TABLE]

where $\alpha$ came from the assumption $\langle u_{y}\rangle/h_{0}=\alpha p(z)/(E_{Y}\delta)$ for the thick top wall and has to be determined by fitting to experiments. Here $\langle u_{y}\rangle$ is the average interface deflection at each fixed-$z$ cross-section. Even though this model has been employed in a lot of later works [25, 24], the unknown fitting parameter $\alpha$ is one of the biggest drawbacks. More recent work has focused on eliminating the fitting parameter, specifically for thinner top walls, plate theory [28, 29] or engineering pressure–displacement models [46] can be used [28, 29] to obtain the hydrodynamic resistance in the deformed microchannel.

Importantly, in the present study, our emphasis is on the thick top wall case, therefore our scaling is different from [28, 29, 26], wherein $\lambda=12(w/t)^{3}\mathcal{P}_{0}/(\overline{E}_{Y}\delta)$ was used. A prefactor $\propto(w/t)^{3}$ shows up in $\lambda$ when the top wall is plate-like, i.e., $t\lesssim w$, which is a consequence of the assumed bending-dominated regime (and results in a very large value of $\lambda$ for thick top walls). Nevertheless, these theories are self-consistent in that the coefficients in the flow rate–pressure drop relation (i.e., the counterpart of $S_{1}$, $S_{2}$ and $S_{3}$ in equation (27) above) become much smaller to balance the large values of $\lambda$. However, in the large-thickness case, we have already shown that $\sigma_{xx}\sim 1/{(\gamma\delta)}^{2}$, meaning that the bending moment in the solid is actually negligible. This observation clearly shows that different solid deformation mechanisms are involved during FSI in microchannels with large versus small top wall thicknesses (compared to the width).

Next, we will give a systematic discussion of the predictions of our FSI theory for microchannels with thick top walls. At steady state, the flow rate $q=const.$, and each cross-section will inflate under the local pressure $p(z)$. The increase in area reduces the local fluid velocity, which introduces a non-constant pressure gradient along the flow-wise $z$-direction. As shown in figure 3(a) for $\lambda=0$ (i.e., when the channel is rigid), the pressure decreases linearly from the inlet to the outlet, and $\mathrm{d}P/\mathrm{d}Z$ is a constant in this case. However, as the top wall becomes softer, with the increase of $\lambda$, the $P(Z)$ profile deviates further from the linear profile, and $\mathrm{d}P/\mathrm{d}Z$ is a decreasing function of $Z$. Accordingly, based on equation (21), the maximum deflection at the fluid–solid interface, $H(0,Z)$, is expected to be concave, as illustrated in figures 3(b) and 4(b).

Equation (27) is applicable in both the steady-state flow-rate-controlled and pressure-drop-controlled flows. As shown in figure 4(a), with controlled flow rate, the total pressure drop is a linear function of flow rate for a rigid channel but a nonlinear function for a soft channel. Furthermore, the pressure drop decreases as the compliance of the channel increases because, under a fixed flow rate, the softer channel will deform more to reduce the flow velocity, and therefore, the pressure losses due to viscosity at each cross-section. In turn, for a pressure-drop controlled flow, the softer channel will allow a higher flow rate, as well as a larger deflection. As shown in figure 4(b) the maximum deflection at the fluid–solid interface also increases with $\lambda$, for a given pressure drop.

Next we compare our theory with previous experimental studies from the literature, namely [21, 26]. Apart from the model (29) proposed in [21], Gervais et al. performed experiment with microchannels with two different Young’s moduli and two different widths. The important parameters are summarized in table 1. Note the thickness for experiment was reported to be larger than 6 mm, but it was numerically shown that 2 mm was thick enough for a sufficiently accurate comparison. Moreover, the undeformed height, $h_{0}$, for the case GEGJ 4 is corrected to 30 $\mu$m instead of the reported 26 $\mu$m based on the value of $\alpha$. We compute $1/(\gamma\delta)^{2}$ in the last column of table 1 and show that our theory is applicable to all the four cases because $1/(\gamma\delta)^{2}\ll 1$ for all data sets. It is important to note that previous theory of microchannel FSI [28, 29] is not applicable, even as an approximation, to any of these cases.

In figure 5, the flow rate–pressure drop relation curves predicted by equation (28) are shown to agree favorably with the experiments. The corresponding predicted maximum displacement at the interface is shown in figure 6. Although some deviations are observed, it is not easy for us to provide a definite reason as to why, due to the lack of information about experimental sources of error in [21]. In the cases GEGJ 2 and 4, there exists an almost constant shift from the experiment, which could be systematic error. For the cases GEGJ 1 and GEGJ 3, the pressure drops predicted by the theory at higher flow rates are larger than the experiments, which would indicate that the theory underestimates the channel deformation at higher flow rates. It is also relevant to note that the worst agreement in figure 5 is for case GEGJ 3, which exhibits the largest deformation in figure 6. Overall, the largest source of uncertainty, however, is the measurement of the undeformed channel height $h_{0}$. Indeed, one reason we have not included the maximum deformation data from [21] in figure 6 is that the error bars are too large to make a meaningful comparison. As we show in the next comparison with the experiments from [26], a small uncertainty in $h_{0}$ can lead to a large effect on the predicted pressure drop.

A further, quantitative, comparison between our theory and the fitting model (29) can be achieved by computing the values of the statistical coefficient of determination $R^{2}$ via least squares [47] for each model, as shown in table 2. Unsurprisingly, the $R^{2}$ values of the model (29) are closer to 1 than those of the present theory because it is a one-parameter best-fit of the experimental data. Nevertheless, the present theory, without any fitting parameters, also give values of $R^{2}\approx 1$, which means that the present fitting-parameter-free theory can capture the physics of the problem as accurately as a fitting model.

More recently, Raj et al. studied the hydrodynamic resistance in microchannels by varying the top wall thickness as well as the Young’s modulus [26]. Six sets of experimental data were reported, with parameters summarized in table 3. A model based on the thick plate assumption was also proposed in [26]. Unfortunately, we have found that the model cannot explain the FSI because the top wall within the present thickness range cannot be regarded as the thick plate (recall section 3). In figure 7, we compare the flow rate–pressure drop relation from our theory to the experiments; once again, favorable agreement is observed. Note that equation (20) does not involve the thickness $t$, which is why only one prediction curve is obtained for all the three thicknesses used in these experiments. The shaded region represents the 5 $\mu$m uncertainty in the underformed channel height reported in [26].

It is well known that microfluidic measurements are highly sensitive to the channel height because $\Delta p\sim\mu lq/(wh_{0}^{3})$ in the lubrication limit [6]. Almost all the experimental data fall into the shaded region, showing that the present theory is able to give quantitative prediction of the hydrodynamic resistance, but perhaps the experiments in [26] were not accurate enough to achieve their goal of addressing the effect of $t/w$. There are slight deviations in the cases RDC 3 and RDC 6 at large flow rates, which might be a consequence of the large elastic deformation in those cases. Also, these two cases have $1/(\gamma\delta)^{2}=0.49$, thus they are already the least favorable ones from the point of view of the limits of applicability of the proposed theory, which requires $1/(\gamma\delta)^{2}\ll 1$. Overall, the agreement between the theoretical predictions and the experiments is quite satisfactory.

We also compare the maximum deflection of the fluid–solid interface at $z=7.5$ mm, with $E_{Y}=0.157$ MPa, as a function of the flow rate. Figure 8 shows a comparison between the theoretical prediction and the experimental data from [26]. The agreement is best for the smaller flow rates. At larger flow rates, the experiment suggest that the deformation saturates, i.e., stops increasing, unlike the theoretical prediction. In this case, we believe that nonlinear elastic effects, which our linear elastic theory cannot capture, begin to dominate. Another two sets of experimental data for $z=15$ mm and $z=22.5$ mm were also provided in [26]. However, we believe that there are potentially some misprints in [26] because according to the previous discussions, the maximum displacement at the fluid–solid interface is expected to be a concave, instead of a convex, function of $z$ (see figures 3(b) and 4(b)). Furthermore, in spite of the deviations of the interface deflection at large $q$, we find that the prediction of the hydrodynamic resistance is still good (see figure 7), which means that the flow rate–pressure drop relation is not sensitive to discrepancies in the maximum deformation, and so it can be pushed to a larger range of flow rates than one might a priori expect.

6 Conclusion

In this work, we presented a theoretical study of the steady-state fluid–structure interaction between a Newtonian fluid and a long and shallow microchannel with a thick compliant top wall. This problem arises in a variety of lab-on-a-chip applications in which soft polymers, such as PDMS, are used to manufacture flow conduits. However, a complete mathematical theory of the “bulging effect” (flow-induced deformation) observed in experiments has remained elusive. To this end, under the lubrication approximation, we appealed to the standard result that the axial velocity profile is parabolic in any cross-section, even if the cross-section varies in the flow-wise direction. Then, importantly, we coupled the fluid and solid mechanics problems through the deformation of the fluid–solid interface (bulging of the microchannel’s top wall).

A scaling analysis of the elastostatics equations (8) for the solid showed that the stresses in the cross-sections perpendicular to the flow-wise direction are dominant and, thus, the 3D solid mechanics problem is simplified to a 2D plane strain problem. Assuming small strains and using the linear theory of elasticity, we show that the top wall’s thickness plays an important role in the stress distribution in solid and, accordingly, has a significant influence on the boundary conditions to be imposed at two lateral surfaces of the top wall. By requiring that width $w$ and thickness $t$ are such that $1/(\gamma\delta)^{2}=(w/t)^{2}\ll 1$, also defined as the large-thickness case in the present study, the top wall deformation was decoupled in the flow-wise direction, allowing us to treat it as a simply supported rectangle at each cross-section. This analysis yielded a self-similar deflection curve at the fluid–solid interface, when scaling the deformation by the pressure. Furthermore, the present analysis showed that the characteristic scale for the interface deformation for the thick-wall problem is independent of thickness as in [21] (but different from the plate-like problem [28, 29, 32]), specifically $\mathcal{P}_{0}w/\overline{E}_{Y}=\mathcal{P}_{0}w(1-\nu^{2})/{E_{Y}}$, which is the expression from [21], $\mathcal{P}_{0}w/E_{Y}$, corrected for a plane strain configuration.

Integrating the flow velocity at a cross-section, we obtained flow rate–pressure drop relation, which deviates from the Poiseuille’s law because it nonlinearly depends on the compliance of the top wall. The results predicted by the present theory agree favorably with the previous experimental studies [21, 26]. While previous theoretical analyses [28, 29, 48, 40] have successfully addressed this type of fluid–structure interaction for thinner plate-like top walls with $t/w\lesssim 1$, the present theory is the first to quantify the hydrodynamic resistance in shallow compliant microchannels with thick top walls such that $(t/w)^{2}\gg 1$.

The present theory is not only fitting-parameter-free but also directly solves for the fluid–solid interface deflection profile without assuming any specific shape. Our theory uncovers the physics hidden in the fitting parameter, $\alpha$, of the widely used model (29), of which many ad-hoc variations have been proposed [25, 39, 26]. The present analysis also provides a clear answer for why the previous plate-theory-based models [28, 29, 41] cannot be pushed to large-thickness regime (even qualitatively) by showing that the bending effects are trivial in the present model. The differences between these theories are also reflected by the different parameter dependencies of the dimensionless numbers quantifying compliance.

Triggered by the observed deviations from the experiments, a remaining question is nonlinear effects, which arise at larger flow rates. One open question is whether large deformations, which occur at higher flow rates, invalidate our assumption of plane strain, perhaps due to significant stretching along the flow-wise direction (see also the discussion of the “RS” data set in [29]), which means that the cross-sectional deformation profiles are no longer decoupled in the flow-wise direction. Additionally, PDMS is known to be porous, thus uptake solvents and swell [49], which might necessitate considering large-deformation poroelastic effects, as in [50, 51].

Another interesting case is when the thickness of top wall of the microchannel is reduced until it can be considered to be “thin.” We provide a preliminary analysis in Appendix, and we argue that the present solid mechanics problem herein should be reduced to a simply supported beam with tension instead of the clamped beam, as in the previous studies using plate theories [28, 29]. Unfortunately, no experimental investigations on this configuration are available to validate the theory against. It would be of interest to conduct direct numerical simulations in future work to further delve into this problem.

\dataccess

The experimental data used is available in the references cited. All plots herein were generated from the equations in the text. No additional data was created in the course of this study. \aucontributeBoth authors contributed to all aspects of this work. \competingWe declare that we have no competing interests. \fundingThis research was supported, in part, by the U.S. National Science Foundation under grant No. CBET-1705637 (to ICC) and a Ross Fellowship from The Graduate School at Purdue University (to XW).

Appendix A Moderate thicknesses and effect on the boundary conditions

Consider a configuration similar to that shown in figure 1 but the thickness of the top wall is not as large. (The side walls are still considered very thick, specifically infinite for the purposes of this discussion.) In this moderate-thickness case, the plane strain assumption is valid as long as the deformation is small enough, but the boundary conditions imposed at the sidewalls are not clear because $\sigma_{xx}$ is not negligible. Specifically, it is not necessarily correct to impose the simply supported boundary conditions from section 3. At the same time, there is no good reason to impose clamped boundary conditions, as the previous studies on microchannels with thinner, plate-like walls [16, 28, 29], because for the geometry considered herein the side surfaces are allowed to deform, while they were assumed to be rigid in those previous studies.

To understand the type of support at the side walls in the moderate-thickness case, we start from the free body diagram in figure 9(b). Then, the reaction forces at the sidewalls are

[TABLE]

where $T_{s}$, $N_{s}$, and $M_{s}$ denote the tension, shear force and the moment respectively. Here $d$ is introduced to represent the point of the reaction force at the bottom of the side solid. Since $h_{0}\ll w$, we expect that $d\ll w$ due to stress concentration. Within the top wall, the resultant tension, $T$, shear force, $N$, and moment, $M$, are expected to be: $T\sim T_{s}$, $N\sim N_{s}$, $M\sim p(z)w^{2}/2+M_{s}\sim p(z)w^{2}(1+\delta^{2}-d/w)/2\sim p(z)w^{2}/2$. Hence, we neglect $M_{s}$ in the following analysis.

Thus, we can consider the configuration from figure 9(c), a slender rectangle subject to pressure at the bottom, and shear and tension forces at the sidewalls. Note that the Airy stress function is still applicable in this case, but it is challenging to solve the corresponding biharmonic equation (9) with the inclusion of tension. Fortunately, the thinness of the structure makes Saint-Venant’s principle applicable, which states that “a local force system has negligible effect on the stress distribution at distances that are large compared with the dimension of the surface where the forces are applied” [52]. Accordingly, the displacement field can be estimated based on a classic engineering model, without knowing the exact details of the stress distribution in solid. Thus, we regard the top wall as a simply supported beam with tension and solve for the displacement field by extending Timoshenko’s beam theory.

Mechanical equilibrium requires that

[TABLE]

where $u^{0}_{y}$ now represents the deflection of the mid-plane of the beam ($y=t/2$). The deformation at the fluid–solid interface is believed to be very close to that of the mid-plane due to the slenderness of the top wall. The corresponding constitutive relations are

[TABLE]

where $I=t^{3}/12$ is the second area moment of the beam cross-section, $\varphi$ represents the rotation of the normal of the cross-section, and $\varkappa$ is the shear correction factor [53]. As before, assuming zero displacement, as well as negligible moment at $x=\pm w/2$, the boundary conditions are

[TABLE]

Equation (31)1 shows that the tension is constant in the cross-section, i.e., $T=T_{s}=p(z)h_{0}$. Then, equations  (31) – (32) can be rewritten in terms of $\varphi$ and made dimensionless:

[TABLE]

Here the constant $\zeta=Tw^{2}/(\overline{E}_{Y}I)=\mathcal{P}_{0}h_{0}w^{2}/(\overline{E}_{Y}I)$ has been introduced to quantify the tension effect. In equation (34), the first term represents the bending effect. The terms in the bracket represents the influence of tension, and the thickness effect is captured by the second term. Given the typical range of parameters for a microchannel, we conclude that the tension cannot be neglected here. Therefore, for the small thickness case, the top wall can no longer be regarded as a simply supported rectangle but, rather, it behaves like a beam with an immovable edge, i.e., simple support plus tension [54].

Soving equations  (31) – (32), the vertical displacement of the mid-plane, which is also approximately the vertical displacement at the fluid–solid interface, is found to be

[TABLE]

where

[TABLE]

If $(t/w)^{2}\ll 1$, then $\mathfrak{u}^{2}(z)\approx\zeta P(Z)/4=p(z)h_{0}w^{2}/(4\overline{E}_{Y}I)$ and the second term in the bracket in equation (35) also vanishes, then equation (35) is reduced to the Euler–Bernoulli beam with tension [54]. However, unlike equations (19)2 and (20) of the large-thickness case, the deflection profile in equation (35) no longer displays self-similarity along the flow-wise direction because $p(z)$ cannot be factored out.

To get a sense of the tension effect, we compare this proposed moderate-thickness theory with other models in figure 10. Note that we compare $\delta H(X,Z)/P(Z)$, with $\delta H(X,Z)=u_{y}^{0}(x,z)/w$, to show the magnitude of the deformation. For the applicability of linear elasticity, we expect that $\delta H(X,Z)\ll 1$. We can see that tension suppresses deformation, compared to the case of simply supported beam. However, the tension is not as restrictive as the clamping considered in [29].

As before, the flow rate–pressure drop relation is obtained by integrating the axial fluid velocity across the cross-section. Then, we rewrite $p(z)$ as a function of $\mathfrak{u}(z)$, so that equation (24) can be written entirely in terms of $\mathfrak{u}$:

[TABLE]

Using separation of variables, the solution of the last ODE can be expressed as a quadrature:

[TABLE]

where $\mathfrak{u}^{\prime}$ is a “dummy” integration variable. Unfortunately, this integral can only be evaluated numerically.