Start-up flow in shallow deformable microchannels

TL;DR

This paper extends the steady flow theory in shallow deformable microchannels to include start-up flow from rest, developing an unsteady lubrication model validated by numerical simulations, revealing insights into transient microfluidic behavior.

Contribution

It introduces a new unsteady lubrication theory for start-up flow in deformable microchannels, considering elastic wall deformation and fluid-solid inertia effects.

Findings

The model accurately predicts transient flow behavior under various conditions.

Self-similar solutions are derived for pressure-controlled start-up flow.

Good agreement with numerical simulations validates the theoretical approach.

Abstract

Microfluidic systems are usually fabricated with soft materials that deform due to the fluid stresses. Recent experimental and theoretical studies on the steady flow in shallow deformable microchannels have shown that the flow rate is a nonlinear function of the pressure drop due to the deformation of the upper soft wall. Here, we extend the steady theory of Christov et al. (2018) by considering the start-up flow from rest, both in pressure-controlled and in flow-rate-controlled configurations. The characteristic scales and relevant parameters governing the transient flow are first identified, followed by the development of an unsteady lubrication theory assuming that the inertia of the fluid is negligible, and that the upper wall can be modeled as an elastic plate under pure bending satisfying the Kirchhoff-Love equation. The model is governed by two non-geometrical dimensionless numbers: a compliance parameter $\beta$, which compares the characteristic displacement of the upper wall with the undeformed channel height, and a parameter $\gamma$ that compares the inertia of the solid with its flexural rigidity. In the limit of negligible solid inertia, $\gamma \to 0$, a quasi-steady model is developed, whereby the fluid pressure satisfies a nonlinear diffusion equation, with $\beta$ as the only parameter, which admits a self-similar solution under pressure-controlled conditions. This simplified lubrication description is validated with coupled three-dimensional numerical simulations of the Navier equations for the elastic solid and the Navier-Stokes equations for the fluid. The agreement is very good when the hypotheses behind the model are satisfied. Unexpectedly, we find fair agreement even in cases where the solid and liquid inertia cannot be neglected.

Tables

Table 1. Table 1: Characteristic scales used in the derivation of the unsteady lubrication model of § 4.3 (upper row), and the dimensionless parameters that govern the elasto-hydrodynamic problem, deduced in § 4.2 (bottom row).

$y$, $x$, $z$	$v_y$, $v_x$, $v_z$ (20a)	Equation (20d)	Equation (20e)	Equation (20f)
$h_0\ll w\ll \ell$, $d\ll w$	$\frac{h_0}{\ell }{V}_{\ell }\ll \frac{w}{\ell }{V}_{\ell }\ll V_{\ell }$	$V_{\ell }\sim \frac{h_0^2{p}_c}{12\mu \ell }$	$u_c\sim \frac{w^4{p}_c}{B}$	$t_c\sim \frac{\mu w^4{\ell }^2}{B{h}_0^3}$

Table 2. Table 2: Physical parameters corresponding to the S4 experiment of  Ozsun et al. ( 2013 ) for two different working liquids, namely water (II), which is the least favorable configuration, and a silicon oil of 500 cSt (I), for which the assumptions of § 4.5 are satisfied. The values of q 0 subscript 𝑞 0 q_{0} and Δ ​ p Δ 𝑝 \Delta p ensure that β ~ = 0.4 ~ 𝛽 0.4 \tilde{\beta}=0.4 , both in the pressure-controlled case and in the flow-rate-controlled situation.

	$h_0$	$w$	$\ell$	$d$	$B$	$\nu$	${\rho }_s$	$\Delta p$	$\tilde{\beta }$
	(mm)	(mm)	(mm)	(mm)	($\mu$J)		(kg m-3)	(kPa)
S4	0.244	1.7	15.5	0.2	$\approx$ 1.6	0.499	970	4.487	0.4

Equations

$$y=h(x,z,t)=h_0+u_y(x,z,t),$$

$$q(z,t)={\int }_{-w/2}^{w/2}{\int }_0^{h(x,z,t)}{v}_z(x,y,z,t)\text{d}y\text{d}x.$$

$$\nabla \cdot \mathbf{v}=0,\text{and}\rho ({\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v})=\nabla \cdot \text{mathsfbi}T,$$

$$\mathbf{v}=0\text{at}y=0,\text{and}x=\pm w/2,$$

$$p=0\text{and}\text{mathsfbi}T\cdot {\mathbf{e}}_z=0\text{at}z=\ell .$$

$$p=\Delta p\text{and}\text{mathsfbi}T\cdot {\mathbf{e}}_z=0\text{at}z=0\text{(pressure-controlled)},$$

$$v_z=v_0(x,y)\text{and}{v}_x=v_y=0\text{at}z=0\text{(flow-rate-controlled)}.$$

$$\mathbf{v}={\partial }_t\mathbf{u}\text{at}y=h(x,z,t).$$

$${\rho }_s{\partial }_t^2\mathbf{u}=\nabla \cdot \mathbf{\sigma },$$

$$\mathbf{\epsilon }=\frac{1}{2}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\text{T}}\right),$$

$$\mathbf{\epsilon }=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^{\text{T}}+\nabla \mathbf{u}\cdot \nabla {\mathbf{u}}^{\text{T}}),$$

$$\mathbf{u}=0\text{at}z=0,z=\ell ,\text{and}x=\pm w/2.$$

$$(\mathbf{\sigma }+\text{mathsfbi}T)\cdot \mathbf{n}=0\text{at}y=h(x,z,t),$$

$$\mathbf{\sigma }\cdot {\mathbf{n}}_{\text{ext}}=0\text{at}y=h(x,z,t)+d(x,z,t),$$

$${\rho }_{s}d{\partial }_t^2{u}_y+B{\nabla }_{xz}^4{u}_y=p\text{at}y=h(x,z,t),$$

$$u_y=0,\text{and}{\partial }_z{u}_y=0\text{at}z=0,\text{and}z=\ell ,$$

$$u_y=0,\text{and}{\partial }_x{u}_y=0\text{at}x=\pm w/2,$$

$${\int }_{-w/2}^{w/2}({\partial }_t{u}_y)\mathrm{d}x+{\partial }_{z}q=0,$$

$$V_{\ell }\sim \frac{h_0^2{p}_c}{12\mu \ell },u_c\sim \frac{w^4{p}_c}{B},\text{and}{t}_c\sim \frac{u_c\ell }{h_0{V}_{\ell }}\sim \frac{12\mu {\ell }^2{w}^4}{B{h}_0^3},$$

$$ϵ=\frac{h_0}{\ell },\delta =\frac{h_0}{w},\alpha =\frac{d}{w},\tilde{\beta }=\frac{w^4{p}_c}{240h_{0}B},ϵRe=\frac{\rho h_0^4{p}_c}{12{\mu }^2{\ell }^2},\tilde{\gamma }=\frac{60^2{\rho }_{s}d{h}_0^{6}B}{{\mu }^2{\ell }^2{w}^4},$$

$$X=\frac{x}{w},Y=\frac{y}{h_0},Z=\frac{z}{\ell },T=\frac{60h_0^{3}B}{\mu w^4{\ell }^2}t,H=\frac{h}{h_0},U_Y=\frac{B}{w^4{p}_c}{u}_y,$$

$$V_X=\frac{\delta }{ϵ}\frac{12\mu \ell }{h_0^2{p}_c}{v}_x,V_Y=\frac{1}{ϵ}\frac{12\mu \ell }{h_0^2{p}_c}{v}_y,V_Z=\frac{12\mu \ell }{h_0^2{p}_c}{v}_z,P=\frac{p}{p_c},Q=\frac{12\mu \ell }{w{h}_0^3{p}_c}q,$$

$$\nabla \cdot \mathbf{V}=0,$$

$$12\left[\frac{{ϵ}^2}{{\delta }^2}\left(\frac{ϵRe}{4\tilde{\beta }}{\partial }_T{V}_X+\frac{ϵRe}{12}\mathbf{V}\cdot \nabla V_X\right)+{\partial }_{X}P\right]={ϵ}^2{\partial }_X^2{V}_X+\frac{{ϵ}^2}{{\delta }^2}{\partial }_Y^2{V}_X+\frac{{ϵ}^4}{{\delta }^2}{\partial }_Z^2{V}_X,$$

$$12\left[{ϵ}^2\left(\frac{ϵRe}{4\tilde{\beta }}{\partial }_T{V}_Y+\frac{ϵRe}{12}\mathbf{V}\cdot \nabla V_Y\right)+{\partial }_{Y}P\right]={ϵ}^2{\delta }^2{\partial }_X^2{V}_Y+{ϵ}^2{\partial }_Y^2{V}_Y+{ϵ}^4{\partial }_Z^2{V}_Y,$$

$$12\left(\frac{ϵRe}{4\tilde{\beta }}{\partial }_T{V}_Z+\frac{ϵRe}{12}\mathbf{V}\cdot \nabla V_Z+{\partial }_{Z}P\right)={\delta }^2{\partial }_X^2{V}_Z+{\partial }_Y^2{V}_Z+{ϵ}^2{\partial }_Z^2{V}_Z,$$

$$\tilde{\gamma }{\partial }_T^2{U}_Y+{\partial }_X^4{U}_Y+\frac{{ϵ}^2}{{\delta }^2}2{\partial }_X^2{\partial }_Z^2{U}_Y+\frac{{ϵ}^4}{{\delta }^4}{\partial }_Z^4{U}_Y=P,$$

$$720{\int }_{-1/2}^{1/2}({\partial }_T{U}_Y)\mathrm{d}X+{\partial }_{Z}Q=0,$$

$$Q(Z,T)={\int }_{-1/2}^{1/2}{\int }_0^{1+\beta U_Y(X,Z,T)}{V}_Z\mathrm{d}Y\mathrm{d}X,$$

$${\partial }_{Z}P={\partial }_Y^2{V}_Z,$$

$$\tilde{\gamma }{\partial }_T^2{U}_Y+{\partial }_X^4{U}_Y=P,$$

$$720{\int }_{-1/2}^{1/2}({\partial }_T{U}_Y)\mathrm{d}X+{\partial }_{Z}Q=0,$$

$$Q(Z,T)={\int }_{-1/2}^{1/2}{\int }_0^{1+\beta U_Y(X,Z,T)}{V}_Z\mathrm{d}Y\mathrm{d}X.$$

$$V_Z=0\text{at}Y=0,\text{and}Y=H(X,Z,T)=1+\beta U_Y(X,Z,T),$$

$$U_Y=0,\text{and}{\partial }_X{U}_Y=0\text{at}X=\pm 1/2,$$

$$\begin{array}{ll}Q=1,& \text{flow-rate-controlled}\\ P=1,& \text{pressure-controlled}\end{array}\}\text{at}Z=0,$$

$$P=0\text{at}Z=1.$$

$$V_Z(X,Y,Z)=-6Y[1+\beta U_Y(X,Z,T)-Y]{\partial }_{Z}P,$$

Full text

1 Introduction

Microfluidic devices allow the manipulation of fluids and objects inside channels whose typical dimensions vary from tens to hundreds of microns. These systems have drastically reduced the working space and the time involved in the applications where they are used. In particular, these devices are routinely used to control multiphase flows, e.g. to generate monodisperse bubbles and droplets, to manipulate immersed soft and hard objects, namely particles, capsules, cells or vesicles, or in diverse applications such as sorting, mixing, drug delivery or mass spectroscopy. The development of microfluidics has facilitated a myriad of applications in a large variety of fields, both in scientific and engineering contexts, e.g. biology (El-Ali et al., 2006), pharmacy and medicine (Rodríguez-Rodríguez et al., 2015) or biomedical research (Sackmann et al., 2014). For more detailed information the reader is referred to the reviews of Stone et al. (2004), Tabeling (2005), Squires & Quake (2005), Whitesides (2006), Bruus (2008) and Anna (2016).

The development of soft lithography (Xia & Whitesides, 1998) played a crucial role in reducing the manufacturing cost and time of these microfluidic platforms. Soft lithography implies the use of highly flexible materials like polydimethylsiloxane (PDMS), since it is cheap, biocompatible, has a low curing time, and is transparent, which facilitates experimental measurements. These materials are cured on a harder substrate, typically glass, which does not deform appreciably under the characteristic overpressures achieved within the channel, in contrast with the PDMS walls, which may experience a substantial deformation. Hence, the use of soft materials in microfluidics naturally gives rise to coupled fluid-structure interaction (FSI) problems under flow (Squires & Quake, 2005; Bruus, 2008; Duprat & Stone, 2015). In fact, many applications take advantage of the deformation of the compliant walls, whose passive or active control allows the development of soft actuators in lab-on-a-chip devices, e.g. valves, pumps, self-regulating components or flow rectifiers, and has also facilitated the development of soft robotics (Ilievski et al., 2011; Majidi, 2014; Elbaz & Gat, 2014; Shepherd et al., 2011; Rus & Tolley, 2015; Polygerinos et al., 2017).

Traditionally, FSI problems have been studied in the context of high-Reynolds-number flows relevant to civil, aeronautical and naval engineering, where aeroelastic and hydroelastic couplings play a crucial role (Païdoussis et al., 2010; Bisplinghoff et al., 2013). In the last few decades, there is a growing attention to FSI problems at small scales due to their ubiquity in nature and in many engineering applications beyond microfluidics. Such is the case of elastocapillarity (Bico et al., 2018), peeling processes (Juel et al., 2018), the flow around swimming bacterial colonies (Lauga, 2016), or around vesicles and blood cells in compliant capillaries (Goldsmith & Skalak, 1975; Secomb et al., 2002), where FSI is crucial to understand the underlying physics. A particularly important context where small-scale FSI problems arise is in biological flows (Fung, 1993a, b, c), e.g. in the pulmonary and respiratory systems (Grotberg, 1994, 2001; Grotberg & Jensen, 2004; Heil & Hazel, 2011). Within this context, a large amount of work has been done to understand the dynamics of fluid-filled elastic tubes that may collapse and buckle due to the transmural pressure (Conrad, 1969; Shapiro, 1977; Cancelli & Pedley, 1985; Pedley & Luo, 1998; Heil et al., 2003). Although most of these studies deal with finite Reynolds numbers, there are also some works dealing with Stokes flow and lubrication theory in this configuration (Heil & Pedley, 1995; Heil, 1997). Furthermore, there are still important open problems concerning hemodynamics in vascular networks, involving vascular remodelling, regulation of blood flow, oxygen transport, the fluid pressure distribution or the shear stress exerted at the compliant vessel walls, which are crucial in the development and detection of cardiovascular diseases such as aneurysms or ischemias, or even in tumor angiogenesis (Goldsmith & Skalak, 1975; Pedley, 1980; Taylor & Draney, 2004; Popel & Johnson, 2005; Cassot et al., 2006; Lasheras, 2007; Sforza et al., 2009).

In the present work we focus on the incompressible start-up flow in a shallow microchannel of rectangular cross-section that is filled with a Newtonian liquid, and whose upper soft wall deforms due to the overpressure needed to induce the flow. It is well known that in the low-Reynolds-number flow of liquids inside rigid channels of constant cross section, the flow rate $q$ is proportional to the pressure drop $\Delta p$. The constant of proportionality is usually referred to as the hydraulic resistance, which only depends on the geometry of the cross-section, on the channel length, and on the fluid dynamic viscosity (Happel & Brenner, 2012). However, several authors have shown through experiments, theory, and simulations that the relationship between $q$ and $\Delta p$ is nonlinear when the deformation of the walls induced by the fluid pressure is not negligible (Gervais et al., 2006; Hardy et al., 2009; Seker et al., 2009; Cheung et al., 2012; Ozsun et al., 2013; Raj & Sen, 2016; Raj et al., 2017; Christov et al., 2018). Under steady flow conditions, these authors found that, for a given imposed pressure drop, the flow rate is larger than that associated with the corresponding rigid channel. Indeed, higher throughputs can be achieved in deformable microchannels due to a decrease in the hydraulic resistance induced by the wall compliance.

The steady lubrication theory developed by Christov et al. (2018) assumes that only the upper wall is deformable. Note that, in most microchannels, the bottom wall is rigid, but the lateral ones are made of the same soft material as the upper one, and thus may also deform. Nevertheless, in shallow geometries the deformation of the lateral walls has a negligible effect outside thin lateral elastic boundary layers, as evidenced by the scaling analysis of Gervais et al. (2006). In addition, the thickness of the lateral walls is typically much larger than the thickness of the upper wall, in which case the lateral wall deformation is much smaller than the upper one. Christov et al. (2018) described the upper wall as a plate under pure bending modelled with the linear Kirchhoff–Love theory (Love, 1888). Previous studies (cf. Gervais et al., 2006; Hardy et al., 2009; Raj & Sen, 2016; Raj et al., 2017) considered a global Hookean relation between the fluid pressure and the spanwise average of the upper wall’s vertical displacement, where the displacement profile is a quadratic function of the spanwise coordinate, instead of the quartic profile predicted by the Kirchhoff–Love plate theory. Introducing the hypothesised Hookean relation into the standard $q$-$\Delta p$ function, these authors deduced a model with one fitting parameter that absorbs the geometric and material constants. This model was able to explain the trends observed in the experiments of Gervais et al. (2006) with thick-walled microchannels. However, Christov et al. (2018) showed that this approximation fails in many configurations, e.g. when the thickness of the top wall is smaller than or comparable to the channel’s width. In contrast, the lubrication theory of Christov et al. (2018) does not have any fitting parameter, and depends only on a compliance parameter $\tilde{\beta}$, that arises naturally from the coupling between the steady fluid flow, described with lubrication theory, and the upper wall displacement, described with the Kirchhoff-Love plate theory. The parameter $\tilde{\beta}$ compares the characteristic displacement of the upper wall with its undeformed height or, equivalently, the characteristic overpressure with the flexural rigidity of the upper wall. Christov et al. (2018) also made direct comparisons of their theory with the experimental data of Ozsun et al. (2013), finding good agreement without fitting parameters. More recently, Shidhore & Christov (2018) have extended these ideas deriving a lubrication model for microchannels with a thicker top wall, whose displacement is modelled with the Mindlin theory accounting for shear stresses. They also performed steady three-dimensional (3D) numerical simulations for both thin and thick upper walls, finding good agreement with their lubrication models, but worse agreement with the experimental data of Ozsun et al. (2013). Gervais et al. (2006) also conducted 3D numerical simulations of the steady flow. However, they did not apply clamped boundary conditions for the top wall displacement at the inlet and outlet of the microchannel, and they did not impose the exact continuity of stresses at the fluid-solid interface.

All the studies mentioned in the previous paragraph focused on steady flow. There are also two previous experimental and theoretical works dealing with unsteady flow in deformable microchannels, namely Dendukuri et al. (2007) and Panda et al. (2009). In particular, both studies considered the stop-flow associated with the relaxation of a top wall initially deflected by the fluid stresses, which induces a squeeze flow towards the inlet and the outlet of the microchannel. To derive their lubrication models, these authors neglected the solid and the liquid inertia, and obtained a nonlinear diffusion equation for the vertical displacement of the upper wall. However, the latter equation is markedly different from the one developed herein, since it is based on a Hookean relation between the vertical displacement of the upper wall and the fluid pressure which is valid in the limit of a very thick wall, but not in the case of a thin plate. In particular, Dendukuri et al. (2007) determined the relaxation time of the upper wall through scaling arguments, finding good agreement with their own experiments with one fitting parameter, equivalent to that introduced by Gervais et al. (2006) in the case of steady flow. In the present work, we show that the model developed by Dendukuri et al. (2007), although an important step forward, fails to describe the unsteady flow for most geometries and wall materials. Another important contribution of Dendukuri et al. (2007) was to show that the characteristic time scale only depends on the geometry of the channel and on the fluid and solid properties, but not on the fluid pressure or flow rate, in agreement with the results developed herein. The lubrication model of Dendukuri et al. (2007) was solved numerically by Panda et al. (2009), and compared with stop-flow experiments performed with thick-walled microchannels. The main limitations of the unsteady lubrication theory developed by Dendukuri et al. (2007) are the same as those of the steady lubrication model of Gervais et al. (2006), whose shortcomings were indicated by Christov et al. (2018).

An unsteady analysis similar to the one presented herein was performed by Elbaz & Gat (2014) for a thin cylindrical soft shell conveying a viscous fluid. These authors identified the characteristic time scale of the unsteady flow, which is equivalent to the one deduced in the present work. In addition, they developed an unsteady lubrication theory neglecting liquid and solid inertia, and deduced a diffusion equation for the fluid pressure. However, the latter equation is linear, since only small values of the compliance parameter, corresponding to small deformations, were considered by Elbaz & Gat (2014).

In the present work, we extend the steady theory of Christov et al. (2018) to account for transient flow, and apply the new framework to the canonical problem of start-up flow from rest. From the theoretical point of view, our main motivation is to provide a framework to tackle FSI problems in laminar internal flows dominated by viscous forces. Indeed, although only the start-up flow is analysed herein for brevity, a similar formalism can be developed to study other transient problems such as stop flows or oscillatory flows. Additional motivation comes from a basic question that, to the best of our knowledge, remains unanswered: what is the start-up time of the flow in a deformable microchannel, and how does it depend on the liquid and solid properties? To that end, we first identify the characteristic hydro-elastic scales and the relevant parameters governing the unsteady flow. Then we develop an unsteady lubrication theory accounting for the solid and liquid inertia, assuming that the upper wall is governed by the Kirchhoff–Love theory in the bending-dominated regime. When the liquid and solid inertia are negligible, we derive a nonlinear diffusion equation for the fluid pressure field. To check the model, we perform 3D direct numerical simulations of the Navier and Navier-Stokes equations for the solid and for the liquid respectively, and compare the results with the quasi-steady lubrication model.

The paper is organised as follows: the flow configuration is described in §2. The mathematical formulation is presented in §3 making use of the Navier and Navier-Stokes equations for the elastic upper wall and the incompressible flow, respectively. In §4 we identify the characteristic scales and the dimensionless parameters governing the flow, and we develop an unsteady lubrication theory for the elasto-hydrodynamic problem, which is further simplified in the quasi-steady limit of negligible liquid and solid inertia. The results are presented in §5, including a comparison between the quasi-steady lubrication model and the 3D numerical simulations. Conclusions are drawn in §6.

2 Flow configuration

As sketched in figure 1, we consider the incompressible start-up flow in a channel of length $\ell$, width $w$ and height $h$, where $h\ll w\ll\ell$, initially filled with a Newtonian fluid of density $\rho$ and dynamic viscosity $\mu$. The overpressure needed to convey the fluid deforms the soft walls, which in turn affects the hydraulic resistance of the channel, giving rise to a coupled fluid-structure problem (Gervais et al., 2006; Weibel et al., 2007; Hardy et al., 2009; Ozsun et al., 2013). Assuming that only the upper wall deforms, here we extend the results of Christov et al. (2018) to unsteady flow.

We adopt a Cartesian coordinate system ($x$, $y$, $z$) as shown in figure 1, and use $\bm{v}=(v_{x},v_{y},v_{z})$ to denote the fluid velocity field and $\bm{u}=(u_{x},u_{y},u_{z})$ to denote the displacement field of the upper wall, of thickness $d(x,z,t)$. We also use $u_{y}(x,z,t)$ to denote the vertical displacement of the lower surface of the upper wall, i.e. the fluid–solid boundary, so that its position is given by

[TABLE]

where $h_{0}$ is the undeformed height of the channel. The displacement is induced by the fluid pressure $p(x,y,z,t)$, which is measured with respect to the outer atmospheric pressure. The flow rate $q(z,t)$ in the $z$-direction is given by the cross-sectional integral of the axial velocity $v_{z}$ as

[TABLE]

For times $t<0$, the fluid is at rest with $p=0$ and $\bm{v}=\bm{0}$, and thus the solid remains undeformed, $\bm{u}=\bm{0}$. For $t>0$, a start-up flow takes place, either due to an imposed inlet flow rate $q(z=0,t)=q_{0}$ (flow-rate-controlled situation) or due to an imposed inlet overpressure $p(z=0,t)=\Delta p>0$ (pressure-controlled situation). We assume that the outlet pressure is $p(z=\ell,t)=0$.

To address the start-up flow we have developed two techniques, namely: 1) 3D numerical simulations of the full Navier-Stokes equations for the flow field and the linear Navier equations for the solid deformation field, and 2) a lubrication theory assuming that the upper wall behaves according to the linear Kirchhoff–Love equation for a plate under pure bending (Love, 1888), neglecting changes in its thickness, $d$.

3 Formulation of the problem

Although the main objective of the present study is to develop an unsteady lubrication theory to describe the start-up flow (see §4), to check its validity we have also performed 3D direct numerical simulations of the Navier equations for a linear elastic upper wall of finite thickness, fully coupled to the Navier-Stokes equations for the incompressible flow of a Newtonian liquid. The fluid velocity field is governed by the continuity and momentum equations,

[TABLE]

where $\mathsfbi{T}=-p\mathsfbi{I}+\mu(\bm{\nabla}\bm{v}+\bm{\nabla}\bm{v}^{\text{T}})$ is the fluid stress tensor. The fluid velocity field must satisfy the no-slip conditions at the rigid walls,

[TABLE]

and atmospheric conditions at the outlet,

[TABLE]

At the inlet, we either impose a fixed pressure $\Delta p$ and no shear stress,

[TABLE]

or a given axial flow profile $v_{0}(x,y)$ corresponding to an input flow rate $q=q_{0}$ (see appendix B) and no tangential flow,

[TABLE]

At the contact surface between the liquid and the soft wall, the continuity of velocities must hold,

[TABLE]

The elastic upper wall satisfies the Navier momentum equation,

[TABLE]

where $\rho_{s}$ is the upper wall material density, $\bm{\sigma}=2\mu_{s}\bm{\varepsilon}+\lambda_{s}\text{tr}(\bm{\varepsilon})\mathsfbi{I}$ is the solid stress tensor, $\bm{\varepsilon}$ is the strain tensor, and $\mu_{s}=E/[2(1+\nu)]$ and $\lambda_{s}=E\nu/[(1+\nu)(1-2\nu)]$ are the two Lamé constants, expressed in terms of the Young modulus $E$ and Poisson ratio $\nu$. The main results of §5 have been obtained assuming that the strain is linear in the displacement gradients,

[TABLE]

but we have also used the complete nonlinear expression of the strain tensor,

[TABLE]

to check whether stretching of the wall has a significant effect. This nonlinear stretching may be significant even when the deformation gradient $\bm{\nabla}\bm{u}$ and strain $\bm{\varepsilon}$ are small and the linear stress–strain relationship holds, such as for a thin elastic sheet that is deflected by an amount $u_{y}$ comparable to or larger than its thickness $d$, which can be modelled using the Föppl–von Kármán equations.

The lateral walls of the solid are clamped, so that $\bm{u}$ satisfies

[TABLE]

At the fluid-solid interface the continuity of stresses must be fulfilled,

[TABLE]

where $\bm{n}$ is the unit normal vector to the liquid-solid interface. Finally, we impose a stress-free condition at the upper surface of the top wall,

[TABLE]

where $\bm{n}_{\text{ext}}$ is the unit normal vector to the upper surface of the soft wall. To perform 3D numerical simulations, the complete system of equations (3)–(14) is expressed in weak form and solved with the finite-element software COMSOL Multiphysics employing the Arbitrary Lagrangian-Eulerian (ALE) method. The details of the numerical techniques employed herein are provided in appendix B.

Figure 2($a$) shows a first comparison of the two calculational frameworks employed herein under steady-state conditions, namely 3D numerical simulations and lubrication theory (Christov et al., 2018). In particular, the results show the pressure drop along the channel $\Delta p$ as a function of the flow rate $q$, which is also used as a validation of the numerical method by comparing with the S4 experiment of Ozsun et al. (2013). These experiments were performed in a PDMS microchannel with water as a working liquid, and whose relevant physical parameters are given in table 2. Our 3D numerical simulations (blue triangles) are in excellent agreement with the experiments (circles), whereas the lubrication approximation (blue solid line) properly captures the nonlinear trend, but slightly overestimates the flow rate, as already shown by Christov et al. (2018). The numerical simulations of Shidhore & Christov (2018) are also shown (red squares), although they were carried out taking a liquid viscosity of $\mu=0.91$ mPa s, which corresponds to a room temperature of approximately 24 oC, whereas ours were computed for $\mu=1$ mPa s, corresponding to a room temperature of 20 oC. The prediction according to lubrication theory for $\mu=0.91$ mPa s is also shown (red solid line). Hence, it is possible to infer that the experiments of Ozsun et al. (2013) took place at a room temperature of approximately 20 oC. Furthermore, our numerical simulations were computed taking $\rho=10^{3}$ kg m*-3*, and $B\simeq 1.6$ $\mu$J, where $B=Ed^{3}/[12(1-\nu^{2})]$ is the bending stiffness of the upper wall. As pointed out by Christov et al. (2018), there is uncertainty in the measurements of the bending stiffness $B$, although the curve $q$-$\Delta p$ is less sensitive to changes in $B$ than to changes in the liquid viscosity $\mu$.

Additionally, we have also considered the nonlinear strain (11) in some simulations where $q$ and $\Delta p$ are relatively high (blue crosses). We observe that there is only a small deviation between the linear and nonlinear results, as will be confirmed in §5, which indicates that the strain is small enough that the linear strain–displacement relationship (10) holds.

A typical 3D simulation is presented in figure 2($b$), which shows the pressure field $p$ in the fluid domain and the vertical solid displacement field $u_{y}$ in the solid domain, under the same conditions as the S4 experiment of Ozsun et al. (2013) for an imposed inlet pressure $\Delta p=6$ kPa. These results correspond to the triangle at the largest value of $q$ in figure 2($a$).

4 Lubrication theory for shallow compliant channels

This section is devoted to the derivation of an unsteady lubrication theory for a slender, shallow and deformable microchannel with rectangular cross-section, using as a starting point the complete set of equations (3)–(14) presented in §3.

4.1 Kirchhoff-Love theory for the upper wall deformation

In developing the lubrication model, the Navier equation (9) will be substituted by an appropriate simplified description, based on plate theory, that takes advantage of the upper wall geometry. In particular, if the maximum displacement of the upper wall is small compared to its thickness $d$, which is constant and smaller than the channel’s width $w$, the dynamics of the plate under pure bending can be described with the Kirchhoff–Love equation (Love, 1888; Howell et al., 2009) with clamped boundary conditions at $x=\leavevmode\nobreak\ \pm w/2$. Therefore, the vertical displacement is now independent of $y$, i.e. $u_{y}=u_{y}(x,z)$, and satisfies

[TABLE]

where $\bm{\nabla}^{4}_{xz}$ is the biharmonic operator in the $(x,z)$ plane, and the bending stiffness $B=\leavevmode\nobreak\ Ed^{3}/[12(1-\leavevmode\nobreak\ \nu^{2})]$ of the upper wall is assumed constant.

4.2 Characteristic scales of the unsteady flow

Here we obtain the characteristic scales and the dimensionless parameters governing the unsteady flow in the coupled elasto-hydrodynamic problem, identifying the conditions under which the set of equations (3)–(8) and (15a)–(15c) can be approximated by either an unsteady or a quasi-steady lubrication model. First, we set the dominant balances that govern the problem by taking advantage of the geometry of the microchannel and its top wall, which both are slender and narrow, i.e. $h_{0}\ll w\ll\ell$, and $d\ll w$. These scales imply the hierarchy $\epsilon\ll\delta\ll 1$, and $\alpha\ll 1$, where $\epsilon=h_{0}/\ell$, $\delta=h_{0}/w$ and $\alpha=d/w$ measure the slenderness and the shallowness of the channel, and the narrow geometry of the plate, respectively. The three dominant balances come from, respectively, the standard lubrication force balance in the $z$-direction $\partial_{z}p\sim\mu\partial_{y}^{2}v_{z}$, the narrow-geometry plate balance in (15a) $B\partial_{x}^{4}u_{y}\sim p$, and the flux balance $\partial_{t}u_{y}\sim v_{y}\sim\epsilon V_{\ell}$. The latter can be also deduced from the Reynolds equation

[TABLE]

where we have used the kinematic condition at the liquid-solid interface, $\bm{v}\bm{\cdot}\bm{n}=\bm{n}\leavevmode\nobreak\ \bm{\cdot}\leavevmode\nobreak\ \partial_{t}\bm{u}$, which yields $\partial_{t}u_{y}+v_{x}\partial_{x}u_{y}+v_{z}\partial_{z}u_{y}=v_{y}$, assuming that the displacement is only in the $y$-direction, $\bm{u}=(0,u_{y},0)$.

Using the balances above, the following relations are obtained:

[TABLE]

where $p_{c}$, $V_{\ell}$, $u_{c}$ and $t_{c}$ are the characteristic pressure, axial velocity, top wall displacement, and time scale, respectively. In particular, the characteristic pressure is $p_{c}\sim\Delta p$ in a pressure-controlled configuration or $p_{c}\sim 12\mu\ell q_{0}/(wh_{0}^{3})$ in a flow-rate-controlled situation, where the factor $12$ is the classical lubrication factor. Note that $t_{c}$ does not depend on $\Delta p$ or $q_{0}$, but only on the geometry and on the fluid and solid properties.

The convective acceleration of the fluid can be neglected compared with the viscous force in (3) when $O(\rho\bm{v}\bm{\cdot}\bm{\nabla}\bm{v})/O(\mu\bm{\nabla}^{2}\bm{v})\sim\epsilon Re=(h_{0}/\ell)\rho h_{0}^{3}p_{c}/(12\mu^{2}\ell)\ll 1$, where $Re$ is the Reynolds number based on the characteristic pressure $p_{c}$ and on the undeformed height $h_{0}$. The local acceleration of the fluid in (3) is negligible when the viscous diffusion time, $\rho h_{0}^{2}/\mu$, is much smaller than the characteristic hydro-elastic time $t_{c}$, i.e. when $O(\rho\partial_{t}\bm{v})/O(\mu\bm{\nabla}^{2}\bm{v})\sim 12\epsilon Re/\beta=\rho h_{0}^{5}B/(\mu^{2}\ell^{2}w^{4})\ll 1$, which is the so-called Womersley number, and where $\beta=u_{c}/h_{0}=w^{4}p_{c}/(h_{0}B)$ is the compliance parameter. Furthermore, the inertia of the top wall can be neglected in (15a) when the characteristic time for which the inertia of the solid affects its displacement, $t_{s}$, is much smaller than the characteristic start-up time involving the deflection of the boundary, $t_{c}$, i.e. when $O(\rho_{s}d\,\partial_{t}^{2}u_{y})/O(B\bm{\nabla}^{4}_{xz}u_{y})\sim\leavevmode\nobreak\ \gamma=(t_{s}/t_{c})^{2}=\leavevmode\nobreak\ \rho_{s}dh_{0}^{6}B/(\mu^{2}\ell^{4}w^{4})\ll 1$. Note that, like $t_{c}$, the dimensionless numbers $\epsilon Re/\beta$ and $\gamma$ do not depend on $\Delta p$ or $q_{0}$, but only on the geometry of the channel and on the solid and fluid properties.

As shown below, a more precise estimate provides $t_{c}\sim\mu\ell^{2}w^{4}/(60Bh_{0}^{3})$, which leads to a modified parameter $\tilde{\gamma}=60^{2}\gamma$. It also proves convenient to define a modified compliance parameter $\tilde{\beta}=\beta/240$, which is ten times smaller than the one used by Christov et al. (2018). As will be shown in §4.3, a more accurate estimate of the local and convective inertia is $\epsilon Re/(4\tilde{\beta})$, and $\epsilon Re/12$, respectively. Taking typical values from the experimental data of Ozsun et al. (2013), who used water as the working liquid, $\tilde{\beta}\sim 1$, $\epsilon Re\sim 1$ and $\tilde{\gamma}\sim 1$. However, it is important to note that $\epsilon Re/(4\tilde{\beta})\propto\mu^{-2}$ and $\tilde{\gamma}\propto\mu^{-2}$, so that both dimensionless parameters rapidly become small as the liquid viscosity increases. Therefore, the flow is governed by six dimensionless parameters, namely $\epsilon$, $\delta$, $\alpha$, $\tilde{\beta}$, $\epsilon Re$, and $\tilde{\gamma}$,

[TABLE]

where the only non-geometrical parameters are $\tilde{\beta}$, $Re$ and $\tilde{\gamma}$. Table 1 summarises the analysis on the characteristic scales presented herein, and the dimensionless parameters governing the problem.

Characteristic times equivalent to those deduced above have been obtained previously. For instance, in the same configuration as ours, Dendukuri et al. (2007) deduced similar scalings under the assumption that the spanwise average of the upper wall’s vertical displacement is linearly proportional to the fluid pressure. However, as shown below, such an assumption leads to a free parameter and cannot describe most channel geometries. Furthermore, following the latter procedure, Tabeling (2005) deduced the characteristic time scale of a deformable cylindrical chamber pumping fluid into a much smaller microfluidic channel, which gives rise to the bottleneck effect, where, due to the deformation of the reservoir, the start-up time increases from minutes to hours. In contrast, Elbaz & Gat (2014) obtained the characteristic start-up time without adjustable parameters, following a procedure similar to that developed in the previous paragraph, but in the case of a cylindrical elastic tube conveying a viscous fluid.

4.3 Non-dimensional formulation

We define the following dimensionless variables

[TABLE]

where $p_{c}=\Delta p$, in a pressure-controlled situation, and $p_{c}=12\mu\ell q_{0}/(wh_{0}^{3})$ in a flow-rate-controlled configuration. Introducing these variables into (3), (15a), and (16) provides the following dimensionless equations

[TABLE]

together with the dimensionless version of the boundary conditions, which are omitted here for simplicity. Here $\bm{V}=(V_{X},V_{Y},V_{Z})$ is the dimensionless velocity vector.

The elasto-hydrodynamic timescale $t_{c}$ is the characteristic time for the channel to transition from the undeformed state to the deformed steady state. For a given channel geometry and working fluid, as the rigidity of the wall increases ($B\to\infty$ and hence $\beta\to 0$) the timescale $t_{c}$ decreases, reflecting the fact that less time is required to inflate the channel to the less deformed final steady state ($u_{c}\to 0$). For a perfectly rigid channel ($\beta=0$), the transition occurs instantaneously (in the absence of fluid inertia). Since $t$ is non-dimensionalised using $t_{c}$ (19), in the limit $\beta\to 0$ the non-dimensional equations obtained below remain regular and tend towards a limiting solution, which we investigate in appendix A.

4.4 Leading-order lubrication model

Assuming $\epsilon\ll\delta\ll 1$, $\epsilon Re\ll 1$ and $\tilde{\beta}\gg\epsilon Re$, the $X$- and $Y$-momentum equation (20b) and (20c) yield $\partial_{X}P=\partial_{Y}P=0$, so that the pressure field is only a function of the axial coordinate $Z$ and time $T$, i.e. $P=P(Z,T)$. Therefore, at leading order, the set of equations (20a)–(20g), yields the following system of nonlinear differential equations,

[TABLE]

Note that the leading-order lubrication equations (21a)–(21d) cannot fulfill all of the boundary conditions of the full set of equations. In particular, the no-slip boundary condition (4) is not satisfied at the side walls $X=\pm 1/2$, and from the kinematic condition (8), only the no-slip condition $V_{Z}=0$ remains to be imposed at $Y=H(X,Z,T)$. In addition, the clamped conditions (15b) cannot be imposed on the upper wall. Therefore, the side-wall boundary conditions for $V_{Z}$, and the inlet and outlet boundary conditions for $U_{Y}$, lead to corrections of the order of $O(\delta)\ll 1$ and $O(\epsilon/\delta)\ll 1$ respectively, as discussed by Christov et al. (2018). Hence, equations (21a)–(21d) must be complemented with the remaining boundary and initial conditions,

[TABLE]

In the present work, instead of tackling the complex set of equations (21a)–(22d), we will just consider the case of negligible solid inertia, corresponding to the limit $\tilde{\gamma}\to 0$.

4.5 The limit $\tilde{\gamma}\to 0$

In the limit $\tilde{\gamma}\to 0$, equations (21a)–(22d) provide, at leading order, the quasi-steady description

[TABLE]

where

[TABLE]

The Reynolds equation (21c) yields, using (23b) and (23c), a nonlinear diffusion equation for $P(Z,T)$,

[TABLE]

subject to the boundary and initial conditions

[TABLE]

Note that if the transient is flow-rate-controlled, the dimensionless inlet pressure is a function of time, $P_{0}(T)=P(Z=0,T)$ in order to impose a constant inlet flow rate $Q_{0}=Q(Z=0,T)=1$, whose expression can be obtained from (23c).

4.6 Self-similar solutions

Under pressure-controlled conditions, equation (25) admits an exact self-similar solution $P=P(\zeta)$ in terms of a rescaled coordinate $\zeta=z/\sqrt{T}$ at early times $T\ll 1$ before the diffusion front has reached the end of the channel. The profile $P(\zeta)$ satisfies the ordinary differential equation and boundary conditions

[TABLE]

where primes indicate derivatives with respect to $\zeta$. Note from (23c) that the corresponding flow rate is of the form $Q=\varphi(\zeta)/\sqrt{T}$, and the result $Q_{0}=\varphi(0)/\sqrt{T}$ is in agreement with the early-time behaviour of the numerical results shown in figure 4($a$) below. Figure 3 shows $P$ and $\varphi$ as a function of $\zeta$ for several values of $\tilde{\beta}$ indicated in the legend. As $\tilde{\beta}$ increases (e.g. bending modulus decreases), the displacement of the top wall increases and hence a higher pressure and flow rate are achieved within the channel.

For $\tilde{\beta}\ll 1$, the self-similar equation (27) becomes the classical diffusion equation with a constant input $P=1$, with solution

[TABLE]

Equation (28) fits the curves of figure 3 when $\tilde{\beta}\ll 1$ and also the early-time trend Q_{0}=1/\sqrt{\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}}T} observed in figure 4($a$). It also agrees with the early-time behaviour of the small-$\tilde{\beta}$ asymptotic solution (35a) derived in appendix A.

Moreover, for $\tilde{\beta}\gg 1$, using the approximation $F(\tilde{\beta}P)\approx(250/3003)\tilde{\beta}^{3}P^{3}$, the equations (25) and (27) become (after a rescaling) the well-known equations for a viscous gravity current and their self-similar form, which have been solved by Huppert (1982).

Although a flow-rate-controlled configuration does not admit an early-time self-similar solution for general $\tilde{\beta}$, it does so in the limits $\tilde{\beta}\to 0$ and $\tilde{\beta}\to\infty$. First, when $\tilde{\beta}\to 0$, equation (25) again becomes the diffusion equation but with constant input flux $-\partial_{Z}P=1$, whose self-similar solution is given by

[TABLE]

The requirement that $F(\tilde{\beta}P)\approx 1$ yields the condition $T\ll\tilde{\beta}^{-2}$ for validity. This result also agrees with the small-$\tilde{\beta}$ asymptotic solution (35b) derived in appendix A, and the result P_{0}=2\sqrt{T/\math@atom{\pi}{\mathchoice{\hbox{\displaystyle\pi}}{\hbox{\textstyle\pi}}{\hbox{\scriptstyle\pi}}{\hbox{\scriptscriptstyle\pi}}}} agrees with the early-time behaviour of figure 4($b$) below. In the opposite limit, $\tilde{\beta}\to\infty$, we again obtain the gravity-current equation but with constant input flux, for which the self-similar solution has $Z\sim\tilde{\beta}^{3/5}T^{4/5}$ and $P\sim\tilde{\beta}^{-3/5}T^{1/5}$ (Huppert, 1982), and the conditions for validity $\tilde{\beta}P\gg 1$ and $Z\ll 1$ yield $\tilde{\beta}^{-2}\ll T\ll\tilde{\beta}^{-3/4}$.

5 Results

This section is devoted to presenting the results obtained with the quasi-steady lubrication theory developed in §4.5, and to compare these results with those extracted from the 3D numerical simulations described in §3. The quasi-steady lubrication equation (25) was solved numerically with a standard finite-difference scheme. In particular, we have computed the fluid pressure distribution $P(Z,T)$, the associated flow-rate distribution, $Q(Z,T)$, and the vertical wall displacement, $U_{Y}(X,Z,T)$, for several values of $\tilde{\beta}$. As shown by Christov et al. (2018), higher flow rates are achieved as $\tilde{\beta}$ increases, corresponding to a larger upper wall deformation and a reduced hydraulic resistance. For instance, in a steady pressure-controlled situation, the dimensionless flow rate is higher for a given pressure drop (Christov et al., 2018),

[TABLE]

since the cross-section increases as the top wall deforms. This result can be observed in figure 4($a$,$b$), which shows the dependence of the inlet flow rate in a pressure-controlled configuration, $Q_{0}=Q(Z=0,T)$ in ($a$), and the inlet overpressure $P_{0}=P(Z=0,T)$ in flow-rate-controlled conditions in ($b$), as a function of time $T$ for several values of $\tilde{\beta}$. The insets show how the solution of Christov et al. (2018) (dashed lines) is reached for $T\gtrsim O(1)$, revealing that $t_{c}$ is indeed the proper time scale. In pressure-controlled configurations $Q_{0}$ decreases with time as the initial condition for the fluid pressure diffuses along the channel and $\partial_{Z}P(Z=0,T)$ decreases in magnitude. In flow-rate-controlled conditions, $P_{0}$ increases with time since the magnitude of $\partial_{Z}P$ decreases as the fluid spreads, so that $F(\tilde{\beta}P)$ must increase to keep $Q_{0}=1$.

To illustrate the diffusion of $P(Z,T)$ and $Q(Z,T)$, and the displacement $\beta U_{Y}(X=0,Z,T)$ along the channel, we show in figure 5 (dashed lines) their time evolution for $\tilde{\beta}=0.4$, under pressure-controlled conditions ($a$) and flow-rate-controlled conditions ($b$). Here, $U_{Y}(X=0,Z,T)=U_{m}(Z,T)=P(Z,T)/384$ is the maximum displacement of the top wall according to equation (23b). The results obtained at the largest time indicated in the legend correspond to the steady solution of Christov et al. (2018).

5.1 Start-up time and comparison with previous models

The start-up time at which the steady-state solution of Christov et al. (2018) is reached, denoted by $T_{d}$, is defined here by the condition $|\partial_{T}Q(Z=1,T_{d})|=10^{-4}$. We have ensured that the value of $T_{d}$ is robust against alternative choices of the dependent variable. For instance, $\max_{Z}|\partial_{T}P(Z,T_{d})|=10^{-4}$ provides very similar values. The start-up time is shown in figures 4($c$,$d$) as a function of $\tilde{\beta}$ in a pressure-controlled and in a flow-rate controlled configuration, respectively, computed with the quasi-steady lubrication model (25), with the 3D numerical simulations of §3, and with the small-deformation ($\tilde{\beta}\ll 1$) asymptotic solutions (35) from appendix A. Figures 4($c$,$d$) reveal that the start-up time decreases as the compliance increases. For instance, the value of $T_{d}$ decreases by a factor of 2 when $\tilde{\beta}$ is increased from [math] to $1.6$ (pressure-controlled) or $2$ (flow-rate-controlled). In the limit $\tilde{\beta}\to 0$, the values of $T_{d}$ reach the corresponding rigid-channel asymptotes, $T_{d}\simeq 0.98$ and $T_{d}\simeq 4.2$, when $\tilde{\beta}\lesssim 0.02$ and $\tilde{\beta}\lesssim 0.01$, respectively. Consequently, the asymptotic solutions developed in §A for $\tilde{\beta}\ll 1$ are valid for $\tilde{\beta}\lesssim 0.01$. In the opposite limit, $\tilde{\beta}\gg 1$, we can obtain the scalings $T_{d}\propto\tilde{\beta}^{-3}$ for a pressure-controlled configuration, and $T_{d}\propto\tilde{\beta}^{-3/4}$ for a flow-rate-controlled situation, from the gravity-current solutions deduced in §4.6, or by setting $u_{c}\sim h_{0}$ in the scaling arguments of §4.2.

The results of the 3D numerical simulations shown in figure 4($c$,$d$) were computed for the geometry of the S4 experiment of Ozsun et al. (2013), different values of $\tilde{\beta}$, and two different combinations of $\epsilon Re$ and $\tilde{\gamma}$, corresponding to two different working liquids (see table 2). Case I corresponds to a silicon oil of viscosity 500 cSt as working liquid, where $\epsilon Re/12=1.9\times 10^{-6}$ and $\tilde{\gamma}=2.1\times 10^{-6}$ (triangles), and Case II corresponds to water, which is the least favourable case for this S4 geometry since $\epsilon Re/12=0.46$ and $\tilde{\gamma}=0.49$ (squares), and thus the inertia of the liquid and of the solid may have influence on the flow. However, we have found fair agreement between the quasi-steady lubrication model with negligible solid and liquid inertia and the 3D numerical simulations in both cases. In particular, the agreement improves for increasing values of $\tilde{\beta}$, indicating that the relative importance of the solid and liquid inertia becomes smaller for larger wall displacements. In particular, the local liquid inertia becomes negligible as $\tilde{\beta}\to\infty$, since $\epsilon Re/(4\tilde{\beta})\to 0$. Morever, $\tilde{\gamma}$ can be expressed in terms of $\tilde{\gamma}\propto\rho_{s}dh_{0}^{2}p_{c}\epsilon^{3}/(12\mu^{2}\ell\tilde{\beta})$, so that the solid inertia becomes negligible as $\tilde{\beta}$ increases for a fixed geometry and liquid.

We have also considered the lubrication model derived by Dendukuri et al. (2007), and later used by Panda et al. (2009), in which the spanwise average of the upper wall’s displacement is assumed to be linearly proportional to the fluid pressure, and thereby the pressure and the displacement fields are only functions of time and the longitudinal coordinate. This type of approximation is usually known as the Winkler foundation (Kerr, 1964). To obtain the dimensionless version of their model we take the characteristic displacement as $u_{c}^{D}=wp_{c}/E$, and the characteristic time as $t_{c}^{D}=\mu\ell^{2}w/(12Eh_{0}^{3})$. Hence, the new dimensionless variables are $t=\mu\ell^{2}w/(12Eh_{0}^{3})\mathcal{T}$ and $u_{y}=(wp_{c}/E)\,\mathcal{U}_{Y}$, and the nonlinear diffusion equation for $P(Z,\mathcal{T})$ reads

[TABLE]

subjected to the same boundary and initial conditions as (25), and where $\mathscr{D}=wp_{c}/(h_{0}E)$ is the associated compliance parameter. Note that, when $\tilde{\beta}\ll 1$ and $\mathscr{D}\ll 1$, (25) and (31) coincide at leading order: $\partial_{T}P=\partial^{2}_{Z}P+O(\tilde{\beta})$, or $\partial_{\mathcal{T}}P=\partial^{2}_{Z}P+O(\mathscr{D})$. However, since the scalings for the pressure and the displacement field are different from our lubrication model, the ratios between the different characteristic scales and the two compliance parameters depend on $\nu$ and on the ratio $w/d$, namely $t_{c}/t_{c}^{D}=12/5(1-\nu^{2})(w/d)^{3}$, $u_{c}/u_{c}^{D}=12(1-\nu^{2})(w/d)^{3}$, and $\tilde{\beta}/\mathscr{D}=(1-\nu^{2})/20(d/w)^{3}$. Taking the S4 experiment of Ozsun et al. (2013) with water as the working liquid (see table 2), and considering $\tilde{\beta}=0.4$ as a typical configuration, then $\mathscr{D}=0.017$, and the dimensional steady-state times predicted by each model are $t_{d}=1.22$ ms (present work) and $t_{d}=1.42$ $\mu$s (Dendukuri et al., 2007; Panda et al., 2009), whose ratio is of the same order as $t_{c}/t_{c}^{D}$ for the parameter values of Case I in table 2. Hence, there is a strong quantitative disagreement between the unsteady lubrication model of Dendukuri et al. (2007), with both our lubrication theory, and the 3D numerical simulations. We thus conclude that the model of Dendukuri et al. (2007) fails to predict the transient flow, especially in microchannels where the thickness of the top wall is smaller than, or of the same order as, the channel width. This situation resembles the shortcomings found in previous steady lubrication models (Gervais et al., 2006; Hardy et al., 2009; Cheung et al., 2012; Raj & Sen, 2016; Raj et al., 2017), in that the fitting parameters that appear in these model have been used even for microchannels with thin upper walls, as pointed out by Christov et al. (2018). Indeed, just like these fitting parameters, $\mathscr{D}$ naively absorbs the geometric and material constants.

5.2 Transient comparison between the quasi-steady lubrication theory and the 3D numerical simulations

To perform a more detailed comparison between the results of the quasi-steady lubrication model (25)–(26c) and the 3D simulations, we have chosen the S4 experiment of Ozsun et al. (2013), whose fixed geometrical and physical parameters are reported in table 2. We have also taken $\tilde{\beta}=0.4$ as a typical value of the compliance parameter, and thus $\Delta p$ and $q_{0}$ (or $V_{\ell}$) are chosen to ensure this value in a pressure-controlled and in a flow-rate-controlled situation, respectively. In the experimental results reported by Ozsun et al. (2013) water was used as working liquid, which corresponds to Case II of table 2 and is an unfavourable case that does not fulfill the lubrication assumptions of §4.5. Note that the latter case corresponds to one of the results of the steady 3D simulations shown in figure 2($a$), which is in excellent agreement with the experiments. We have also considered a more favourable case, namely Case I of table 2, so that the hypotheses behind (25)–(26c) are satisfied.

Figures 5 and 6 show the flow rate $Q$ (upper row), the pressure distribution $P$ (middle row) and the displacement of the wall $\beta U_{m}$ (bottom row) as functions of $Z$ at different times indicated in the legend, obtained from the lubrication theory (dashed lines) and from the 3D simulations (solid lines). In both figures, the configuration is pressure controlled in the left column, and flow-rate controlled in the right column. In the 3D simulations, the pressure drop is evaluated along the line $X=0,Y=1/2$, while the solid deformation field is evaluated at the fluid-solid interface and $X=0$. The control parameters correspond to Case I in figure 5 and to Case II in figure 6.

In Case I, the agreement in the time evolution of $Q$, $P$ and $U_{m}$ is fairly good both in the pressure-driven and in the flow-rate-controlled configurations, as evidenced by figure 5. In particular, the self-similar solution $P(\zeta)$, $\varphi(\zeta)$ given by (27) correctly describes the rescaled numerical solution, as shown by the insets in figure 5($a$). In this case, since $\epsilon Re/12,\epsilon Re/(4\tilde{\beta}),\tilde{\gamma}\ll 1$, the largest source of error is probably the fact that the lubrication approximation does not satisfy the no-slip condition for $\bm{v}$ at the lateral walls, nor does it satisfy the clamped condition for $\bm{u}$ at the inlet, which gives rise to an elastic boundary layer where the largest disagreement takes place, especially for $U_{m}$. However, as pointed out by Christov et al. (2018), its influence is confined to a region of length $O(\epsilon/\delta)\ll 1$. Additionally, the lubrication approximation is not able to capture the early-time oscillations experienced by $P$ and $U_{m}$ close to the inlet. These travelling waves are always present even when $\epsilon Re$ and $\tilde{\gamma}$ are exactly zero. Hence, a possible explanation might be that the derivatives $\partial^{4}_{z}u_{y}$ and $\partial^{2}_{z}\partial^{2}_{x}u_{y}$ are significant at early time since its characteristic length scale is $Z\sim 0.05-0.1$, which is of the same order as $\delta$ for the geometry considered in Cases I and II. Therefore, the $Z$-scale is initially small as the fluid spreads within the channel causing the oscillations, which eventually disappear as the fluid propagates and the $Z$-scale increases. An analogous phenomenon has been observed by Lister et al. (2013), where they found travelling-wave solutions for the peeling of an axisymetric elastic sheet.

In Case II, although the values of $\epsilon Re/12$, $\epsilon Re/(4\tilde{\beta})$, and $\tilde{\gamma}$ are not strictly small, and the lubrication hypotheses are not satisfied, the agreement between the 3D simulations and the quasi-steady approximation is better than might be expected, as shown by the results in figure 6. However, in this case, the amplitude and the dissipation time of the early-time oscillations of $P$ and $U_{m}$ are larger, and they also propagate downstream to larger values of $Z$. This behaviour breaks the self-similarity of $P$ and $Q$ in the pressure-controlled configuration.

Finally, to provide a better illustration of the agreement between the quasi-steady lubrication model and the 3D numerical simulations, we have also computed the time-dependent evolution of the longitudinal velocity $V_{Z}$ and the vertical displacement $\beta U_{Y}$ at a longitudinal station close to the inlet, $Z=0.25$, under flow-rate-controlled conditions, and for the values of Case II in table 2, i.e. the least favourable configuration, but for $\tilde{\beta}=4.8$, in order to also test the validity of the linear-strain approximation (10). These conditions correspond to an inlet flow rate of $q_{0}=429$ ml min*-1*, whereas the maximum flow rate reported in the experiments of Ozsun et al. (2013) is $q_{0}=50$ ml min*-1* (see figure 2). Figure 7 displays six different time snapshots, the last one corresponding to the steady state, showing $V_{Z}$ and $\beta U_{Y}$ obtained from the quasi-steady lubrication model (25)–(26c) in the column ($a$–$f$), from the 3D numerical simulations, using the linear strain (10), in the column ($g$–$l$), and using the nonlinear strain, in the column ($m$–$r$). Taking into account that the hypotheses behind the lubrication model (25)–(26c) are not strictly satisfied, the overall agreement between the quasi-steady model and the 3D simulations during the whole start-up transient is quite good. However, again, there are marked differences between both approaches, e.g. in the no-slip condition at the lateral bounding walls, which the lubrication model cannot fulfill, or the constant thickness and the unidirectional displacement of the top wall considered by the bending-dominated Kirchhoff–Love theory, which do not apply in this configuration where the top wall thickness is comparable to the channel width. Furthermore, the quasi-steady lubrication model slightly underestimates the axial velocity of the fluid and the vertical displacement.

Comparing the 3D numerical results ($g$–$l$) and ($m$–$r$) for linear (10) and nonlinear (11) strains, respectively, we find that in this case, even though $\tilde{\beta}$ is relatively large and the deflection is comparable to the thickness of the elastic wall, the nonlinear stretching of the wall has a negligible effect on its elastic response, and the linear-strain model is adequate.

6 Conclusions

In this paper we have studied the start-up flow in a shallow rectangular microchannel with a deformable top wall, considering both pressure-controlled and flow-rate-controlled conditions. To that end, we have developed an unsteady lubrication model, where the top wall is modelled with the Kirchhoff–Love plate theory in the bending-dominated limit. To derive this simplified model we have first identified the characteristic scales and the dimensionless parameters governing the hydro-elastic problem showing, in particular, that the characteristic start-up time only depends on the geometry of the channel and on the solid and fluid properties, but not on the characteristic pressure and flow rate. When the solid and liquid inertia are negligible, the lubrication model is quasi-steady and reduces to a nonlinear diffusion equation for the fluid pressure field, whose only dimensionless parameter is the compliance parameter $\tilde{\beta}$ and that, under pressure-controlled conditions, admits a self-similar solution.

To check the validity of the hypotheses behind the lubrication model in the limit of negligible solid and liquid inertia, we have conducted 3D numerical simulations of the complete Navier and Navier-Stokes equations for the solid and for the fluid, respectively. In particular, as a basis for comparison, we have selected a microchannel whose geometry corresponds to the S4 experiment of Ozsun et al. (2013), and two working liquids. First, we have considered a silicon oil of 500 cSt dynamic viscosity, which fulfills the hypotheses, and we have found excellent agreement between the pressure, displacement, and flow rate predicted by the quasi-steady lubrication model, and those obtained from the 3D numerical simulations. In the second case we have considered water as working liquid, for which the liquid and solid inertia are not negligible, although its influence is moderate. In this case we have also obtained fair agreement. We have also derived a leading-order asymptotic solution in terms of a regular expansion in the compliance parameter $\tilde{\beta}$, which properly captures the transient dynamics of the microchannel when $\tilde{\beta}\ll 1$.

We have also computed the start-up time for several values of $\tilde{\beta}$, comparing the values predicted by our lubrication model with those obtained from the model of Dendukuri et al. (2007) and from the 3D numerical simulations. In the two flow configurations considered herein, we have obtained good agreement between the simulations and our model, but not with the lubrication approximation of Dendukuri et al. (2007). The reason is that these authors assume that the pressure scales as $p_{c}\sim u_{c}E/w$, which is only valid if the elastic wall is large enough, but not when the thickness of the wall is smaller than, or of the same order as the channel width, where the appropriate scaling is $p_{c}\sim u_{c}B/w^{4}$.

There are many extensions of the present work that deserve further effort. First of all, the inertial corrections associated with finite values of $\tilde{\alpha}$ and $\tilde{\gamma}$ in the lubrication equations (21a)–(22d) should be studied. In addition, other unsteady processes like pulsatile flows should be addressed. For intermediate and large wall thicknesses, the bending-dominated Kirchhoff–Love plate theory fails, and thus stretching and shear have to be included in the modelling, which could be addressed along the lines of Shidhore & Christov (2018) for steady flow. Finally, the case of a microchannel embedded in an elastic half-space should also be studied, since it occurs frequently in applications.

Acknowledgements.

Acknowledgements

The authors are grateful to Javier Rivero-Rodríguez and Benoit Scheid for key numerical advice, to Ivan C. Christov for pointing out a mistake in figure 1 of an earlier version of the manuscript, and to Ramón Zaera for helpful discussions. AM-C and AS thank the Spanish MINECO, Subdirección General de Gestión de Ayudas a la Investigación, for its support through projects DPI2014-59292-C3-1-P and DPI2015-71901-REDT, and the Spanish MCIU-Agencia Estatal de Investigación through project DPI2017-88201-C3-3-R. These research projects have been partly financed through FEDER European funds. AM-C also acknowledges support from the Spanish MECD through the grant FPU16/02562 and to its associated program Ayudas a la Movilidad 2018 during his stay at the Complex Fluids Group in Princeton. HAS thanks the NSF for support via CMMI-166-1672 and through Princeton University’s Material Research Science and Engineering Center DMR-1420541.

—————————–

Appendix A The small-compliance limit, $\tilde{\beta}\ll 1$

To study the limit $\tilde{\beta}\ll 1$, we expand the pressure field in $\tilde{\beta}$ as

[TABLE]

At leading order, using $F(\tilde{\beta}P)\approx 1$, the governing equation (25) simplifies to the diffusion equation,

[TABLE]

with boundary and initial conditions

[TABLE]

The solutions can be found using (for example) a Fourier expansion, and are given by

[TABLE]

The first-order correction $P^{(1)}$ satisfies the equation

[TABLE]

with suitable boundary conditions, and can in principle be calculated in the same way, but the leading-order result is sufficient to verify our numerical results for $\tilde{\beta}\ll 1$.

Figure 8 shows a comparison of the inlet flow rate $Q_{0}(T)$ in a pressure-controlled configuration ($P_{0}=1$) and the inlet overpressure $P_{0}(T)$ in a flow-rate-controlled configuration ($Q_{0}=1$), between numerical computations of (25)–(26c) and the corresponding leading-order solutions (35a) and (35b), respectively. These solutions have been used to obtain the start-up times $T_{d}$ in the limit $\tilde{\beta}\ll 1$, shown in figure 4($c$,$d$) (dashed line). Note that the leading-order asymptotic solutions work reasonably well when $\tilde{\beta}\lesssim 0.1$.

Appendix B Numerical implementation

In this appendix we describe the numerical techniques used to implement the system (3)–(14). All the equations are written in weak form by means of the corresponding integral scalar product, defined in terms of test functions for the pressure, velocity, and displacement fields, i.e. $\tilde{p}$, $\tilde{\bm{v}}$, and $\tilde{\bm{u}}$, respectively. By using Green identities we finally obtain an integral bilinear system of equations for the set of variables and their corresponding test functions. Equation (3) reads in weak form

[TABLE]

where $\Omega_{f}$ is the deformable 3D fluid domain, and $\Gamma_{f}^{i}$ are the boundary surfaces with their corresponding normal vectors $\bm{n}_{f}^{i}$. The flux integral is set to zero at the fluid-solid interface, since the continuity of stress is imposed in the weak-form Navier equation. At the lateral and lower walls the no-slip condition for the velocity field (4) is taken into account by imposing $\bm{\tilde{v}}=\bm{0}$ at the corresponding boundaries. At the outlet we impose $\tilde{p}=0$ as boundary condition for the pressure field, and also the stress-free boundary condition (5). Furthermore, at the inlet, in a pressure-controlled situation we impose a non-homogeneous Dirichlet boundary condition for the pressure and a stress-free boundary condition (6). On the other hand, when the microchannel is flow-rate controlled we impose a normal velocity $v_{0}(x,y)$ corresponding to flow in a rigid channel with the desired flow rate, i.e.

[TABLE]

and homogeneous Dirichlet boundary conditions for the tangential velocity (7).

The Navier equation (9) reads in weak form

[TABLE]

where $\Omega_{s}$ is the 3D solid domain, and $\Gamma_{s}^{i}$ the boundaries with $\bm{n}_{s}^{i}$ the corresponding normal vector. The four clamping conditions (12) are imposed as $\tilde{\bm{u}}=\bm{0}$ in the flux integral, whereas the continuity of stress (13) at the contact interface and the stress-free condition (14) at the outer wall are imposed as natural boundary conditions, which read, respectively

[TABLE]

[TABLE]

Moreover, the continuity of velocity (8) is imposed as a weak constraint.

The equations are discretised using Taylor-Hood tetrahedral elements for the pressure and velocity fields, and second-order Lagrange elements for the displacement field, which ensures numerical stability. To account for the deformation of the domain we use the ALE method implemented in COMSOL Multiphysics, where the mesh elements in the solid domain $\Omega_{s}$ move with an imposed displacement given by $\bm{u}$, whereas in the fluid domain $\Omega_{f}$ they move according to the Laplace equation for the change of variable between the material and the spatial frames (Rivero-Rodríguez & Scheid, 2018; Rivero-Rodriguez & Scheid, 2019). As for the time-stepping, we employ a 4th-order variable-step BDF method, or an implicit generalised-alpha method when the inertia of the solid becomes relevant (Case II). The relative tolerance of the nonlinear method is always set to $10^{-5}$.