Capillary-flow dynamics in open rectangular microchannels

Abstract

Spontaneous capillary flow of liquids in narrow spaces plays a key role in a plethora of applications including lab-on-a-chip devices, heat pipes, propellant management devices in spacecrafts and flexible printed electronics manufacturing. In this work we use a combination of theory and experiment to examine capillary-flow dynamics in open rectangular microchannels, which are often found in these applications. Scanning electron microscopy and profilometry are used to highlight the complexity of the free-surface morphology. We develop a self-similar lubrication-theory-based model accounting for this complexity and compare model predictions to those from the widely used modified Lucas–Washburn model, as well as experimental observations over a wide range of channel aspect ratios $\lambda$ and equilibrium contact angles $\theta _0$. We demonstrate that for large $\lambda$ the two model predictions are indistinguishable, whereas for smaller $\lambda$ the lubrication-theory-based model agrees better with experiments. The lubrication-theory-based model is also shown to have better agreement with experiments at smaller $\theta _0$, although as $\theta _0\rightarrow {\rm \pi}/4$ it fails to account for important axial curvature contributions to the free surface and the agreement worsens. Finally, we show that the lubrication-theory-based model also quantitatively predicts the dynamics of fingers that extend ahead of the meniscus. These findings elucidate the limitations of the modified Lucas–Washburn model and demonstrate the importance of accounting for the effects of complex free-surface morphology on capillary-flow dynamics in open rectangular microchannels.

1. Introduction

Capillary flow is the spontaneous wicking of liquid in narrow spaces without the assistance of, or even in opposition to, external forces such as gravity. This phenomenon has been investigated since the early twentieth century and has been exploited for a diverse range of applications including lab-on-a-chip devices (Olanrewaju et al. 2018), heat pipes (Faghri 1995), propellant management devices in spacecrafts (Levine et al. 2015) and fabrication of flexible printed electronics (Cao et al. 2018; Jochem et al. 2018).

Early studies focused on understanding the physical mechanism driving spontaneous capillary flow in capillary tubes. Lucas (1918) and Washburn (1921) appear to have been the first to propose theoretical models describing the meniscus position $\tilde {z}_m$ as a function of time $\tilde {t}$ for flow of a Newtonian liquid in cylindrical capillaries. Lucas (1918) assumed the flow is driven by the capillary pressure gradient caused by the circular-arc meniscus front, while Washburn (1921) also included hydrostatic pressure gradients and an imposed pressure difference between the two ends of the capillary. For a horizontal capillary tube open at both ends, an analytical solution $\tilde {z}_m=\sqrt {\tilde {k} \tilde {t}}$ is obtained, commonly referred to as the Lucas–Washburn relation, where $\tilde {k}$ is known as the mobility parameter and depends on the cylinder radius, liquid viscosity, surface tension and contact angle. The mobility parameter $\tilde {k}$ can be thought of as a diffusion coefficient driving the growth of the liquid interface.

Numerous studies extended the theoretical work of Lucas (1918) and Washburn (1921) by including inertial (Quéré 1997; Rideal 1922; Bosanquet 1923), dynamic contact angle (Siebold et al. 2000; Popescu, Ralston & Sedev 2008; Ouali et al. 2013) and surface roughness (Ouali et al. 2013) effects. Additionally, these theoretical models have been extensively compared to experiments (Rideal 1922; Fisher & Lark 1979; Ichikawa & Satoda 1994; Quéré 1997; Ichikawa, Hosokawa & Maeda 2004; Ouali et al. 2013), confirming the $\tilde {z}_m \sim \tilde {t}^{1/2}$ scaling.

Due to breakthroughs in lithographic fabrication techniques, open microchannels with various cross-sectional geometries can be fabricated easily and inexpensively, including rectangular (Yang et al. 2011; Sowers et al. 2016; Lade et al. 2018; Kolliopoulos et al. 2019), trapezoidal (Chen 2014), U-shaped (Yang et al. 2011) and V-shaped (Rye, Yost & O'Toole 1998; Mann et al. 1995; Rye, Mann & Yost 1996; Yost, Rye & Mann 1997) cross-sections. The lack of a top provides access to the inside of the channel, and has been exploited in applications such as capillary micromoulding and microfluidics. Some studies have generalized the Lucas–Washburn relation to arbitrary cross-sectional geometries (Ouali et al. 2013; Berthier, Gosselin & Berthier 2015). However, predictions of the modified Lucas–Washburn models for open capillaries have resulted in varying agreement with experiments (Yang et al. 2011; Ouali et al. 2013; Chen 2014; Sowers et al. 2016; Kolliopoulos et al. 2019). This is because the mechanism for capillary flow in open channels is more complex than that for closed channels. While for closed channels the force driving the flow is due to the pressure gradient caused by the circular-arc meniscus front, for open channels the additional free surface also contributes to driving the flow (this will be discussed in more detail when presenting figure 4).

The additional contribution of the free-surface curvature to capillary flow has been theoretically and experimentally investigated primarily for V-shaped channels (Mann et al. 1995; Romero & Yost 1996; Rye et al. 1996; Yost et al. 1997; Rye et al. 1998; Weislogel & Lichter 1998; Weislogel 2012). However, while the most widely used open-channel cross-sectional geometry is rectangular (Olanrewaju et al. 2018), previous theoretical studies have only considered capillary flow in open rectangular channels for liquids with contact angles of $\theta _0=0^{\circ }$ and large channel aspect ratios $\lambda =H/W$ (height/width) (Tchikanda, Nilson & Griffiths 2004; Nilson et al. 2006), or reported three-dimensional simulations using the volume-of-fluid method to study the effects of gravity on capillary rise in open rectangular channels (Gurumurthy et al. 2018).

In open rectangular channels the free-surface morphology is more complex than in V-shaped channels. From the channel inlet to the meniscus front the upper meniscus spans the entire channel width. However, at the meniscus front the flow splits into the channel corners provided the equilibrium contact angle $\theta _0<{\rm \pi} /4$ (Concus & Finn 1969). This splitting of the flow leads to filaments or fingers extending ahead of the meniscus front and influencing its propagation. Such a transition is not observed in V-shaped channels.

In this work we use a combination of experiment (§ 2) and theory (§ 3) to study capillary-flow dynamics in open rectangular channels. This is achieved by developing a self-similar lubrication-theory-based model (§ 3.2), and comparing model predictions to the modified Lucas–Washburn (MLW) model (§ 3.1) and complementary flow visualization experiments. We investigate the effects of the complex free-surface morphology on the flow dynamics over a wide range of channel aspect ratios $\lambda$ and equilibrium contact angles $\theta _0$ (§ 4.1) and identify limitations of the MLW model (§ 4.2). Finally, we show good agreement between lubrication-theory-based model predictions of the finger dynamics and experiments (§ 4.3).

2. Capillary-flow experiments

Experiments with a non-volatile liquid are used to study capillary flow in open rectangular microchannels. Flow visualization is used to track the meniscus front and a combination of scanning electron microscopy (SEM) and profilometry is used to characterize the effect of channel aspect ratio on the free-surface morphology.

2.1. Channel fabrication and materials characterization

2.1.1. Fabrication of master pattern

Traditional microfabrication techniques were used to form silicon master patterns of capillary channels. A 10.2 cm diameter silicon wafer was cleaned in an oxygen asher (Technics Oxygen Asher) for 5 min with an oxygen flow of 200 SCCM and RF power of 250 W. MicroChem SU-8 2010 negative tone photoresist was spin-coated onto the wafer at 300 rpm for 5 s and 1000 rpm for 30 s, followed by edge-bead removal with MicroChem EBR PG. These coating conditions target a $20\ \mathrm {\mu }\textrm {m}$ layer thickness. Fabrication of capillary channels using SU-8 was chosen because it gives smother sidewalls, sharper bottom corners and a flatter channel bottom than deep reactive-ion etching. The resist was soft-baked on a hot plate at $95\,^{\circ }\textrm {C}$ for 4 min. The photoresist was exposed through a photomask using a Karl Suss MA6 contact mask aligner in soft contact mode for 12.5 s with a $50\ \mathrm {\mu }\textrm {m}$ gap to define the capillary channels. Measurement gradient marks were included in the master pattern to facilitate tracking of the capillary flow. The wafer was then baked at $95\,^{\circ }\textrm {C}$ for 4 min. The exposed wafer was developed in propylene glycol momomethyl ether acetate (Sigma Aldrich) and rinsed with isopropanol. The resist was then hard baked at $150\,^{\circ }\textrm {C}$ for 30 min and an anti-stick fluorinated monolayer was formed by placing the dried wafer in a reduced-pressure chamber with trichloro(1H,1H,2H,2H-perfluorooctyl)silane (Sigma Aldrich) vapour overnight. The resulting microchannel height was $22.5\ \mathrm {\mu }\textrm {m}$, measured with a KLA Tencor P16 surface profilometer.

2.1.2. Substrate fabrication

Capillary channels were prepared by first casting a silicone stamp (Sylgard 184) over the master pattern, curing the stamp and then using the stamp to imprint UV-curable adhesive (Norland Products NOA68 or NOA73) as explained in Kolliopoulos et al. (2019). Briefly, the UV-curable adhesive was coated on glass slides and then the silicone stamp was pressed into the adhesive. The adhesive was solidified by exposure to 365 nm UV light (Honle UV Spot 100) at $30\ \textrm {mW}\ \textrm {cm}^{-2}$ for 270 s. The stamp was then delaminated from the prepared capillary channels and the channels were inspected with a digital microscope for defects. Any channels with defects were not used for capillary-flow experiments. The microchannel length and height were 30 mm and $22.5\ \mathrm {\mu }\textrm {m}$, respectively. Microchannel widths were 17, 25, 50, 75, 100 and $200\ \mathrm {\mu }\textrm {m}$. The reservoir radius was 3 mm. A schematic of the microchannel geometry and a SEM image of a $100\ \mathrm {\mu }\textrm {m}$ wide and $22.5\ \mathrm {\mu }\textrm {m}$ deep channel are shown in figure 1.

Figure 1. (a) Schematic of microchannel connected to reservoir and (b) SEM image of channel cross-section for $\lambda =H/W=0.45$ (scale bar: $20\ \mathrm {\mu }\textrm {m}$).

2.1.3. Materials characterization

The non-volatile test liquids chosen for capillary-flow experiments included UV-curable adhesive (NOA74, Norland Products), silicone oil (DC-704, Dow Corning Corporation), mineral oil (Sigma-Aldrich) and propylene glycol (Froggy's Fog). Shear viscosity $\mu$ was measured using a stress-controlled rheometer (AR-G2, TA Instruments) with a stainless steel cone-and-plate geometry (40 mm, $2^\circ$ cone angle). Surface tension $\sigma$ was measured using a Krüss DSA-30 digital tensiometer. A Krüss goniometer was used to measure equilibrium contact angles $\theta _0$ on flat test substrates prepared in the same way as the capillary channels. Density values were obtained from the manufacturer specifications. The physical properties and equilibrium contact angles of the test liquids are shown in table 1. Note that all liquids have $\theta _0 < 45^{\circ }$.

Table 1. Physical properties and equilibrium contact angles of test liquids.

$^{a}$Solid NOA73 substrate.

$^{b}$Solid NOA68 substrate.

2.2. Experimental methods

2.2.1. Capillary-flow visualization

The experimental investigation of capillary flow was conducted with the apparatus depicted in figure 2(a). Capillary channels were placed on a custom-built motorized stage assembly which was lit from below through the transparent stage and substrate. A controlled volume of the test liquid was placed into the reservoir attached to the capillary channel using a Nordson EFD ValveMate 7100 drop dispensing system with a 25 GA Nordson EFD tip mounted above the reservoir. Sufficient liquid was deposited at the centre of the reservoir to fully fill the reservoir. After deposition, a programmed microstepping motor (Automation Direct STP-MTRD-23042RE) moved the stage assembly and ensured the liquid front remained in the field of view, allowing for visualization of longer flow distances compared with prior studies (Yang et al. 2011; Sowers et al. 2016; Lade et al. 2018; Kolliopoulos et al. 2019). A high-speed camera (Photron Fastcam-Ultima APX) with a Micro-Nikkor 105 mm lens, Nikon PN-11, Nikon PK-13 and Kenko 20 mm and 36 mm extension tubes and a Kenko N-AFD $2\times$ Teleplus MC7 lens was used to visualize the flow at 60 frames per second. Flow was recorded until the liquid meniscus reached the end of the 30 mm long channel or until the maximum recording time of the camera ($\sim$400 s) was reached. Experiments were conducted at ambient conditions ($23 \pm 1\,^\circ \textrm {C}$). Flow visualization experiments were analysed using ImageJ software. A minimum of four trials were conducted for a given channel aspect ratio and test liquid. The meniscus-position time evolution $\tilde {z}_m(\tilde {t})$ was averaged over all trials and the maximum and minimum $\tilde {z}_m(\tilde {t})$ were used for the range of experimental results. The meniscus-position time evolution $\tilde {z}_m(\tilde {t})$ for different channel aspect ratios $\lambda$ using NOA74 as the test liquid can be seen in figure 2(b). Results using the other test liquids are reported in § 4.2.

Figure 2. (a) Schematic of experimental apparatus. (b) Meniscus position $\tilde {z}_m$ as a function of time $\tilde {t}$ for different channel aspect ratios $\lambda$ with NOA74. The solid lines and shaded areas represent the average and range of experimental results, respectively. The filled symbols indicate the meniscus position and flow time at which samples were cured to obtain SEM and profilometry images (figure 3). For $\lambda =0.45$ the shaded region is so small it cannot be seen.

2.2.2. Free-surface profile characterization

The following experiments were conducted to investigate the effect of channel aspect ratio $\lambda$ on the free-surface morphology. A UV-curable liquid (NOA74, Norland Optical Adhesives) was deposited in the reservoir connected to the microchannel and allowed to flow along the channel length. The flow was terminated at a desired time by exposing the NOA74 to a high-intensity UV light source (Omnicure S1500A with a custom light guide) at a UV intensity of approximately $1.6\ \textrm {W}\ \textrm {cm}^{-2}$. The liquid was fully cured in $<$2 s, but the flow terminates well before full solidification, so the process essentially creates a snapshot of the free-surface profile at a given time. The position of the meniscus front and the time at which curing occurred are represented by the filled symbols in figure 2(b) for two channel aspect ratios $\lambda$.

After solidification, the free-surface profile was measured with a stylus profilometer (KLA-Tencor P16) by making repeated scans across the channel width. The samples were then coated with a conductive gold film and the region near the meniscus front was imaged with a SEM instrument (JEOL JSM-6010PLUS/LA) in secondary electron imaging mode with a sample rotated $40^{\circ }$ about the $\tilde {z}$ axis (figure 3f). The SEM images and profilometry scans for channels of aspect ratio $\lambda =0.45$ and $\lambda =0.225$ corresponding to the filled symbols in figure 2(b) are shown in figure 3.

Figure 3. Profilometry scans and SEM images of cured NOA74 in channels of aspect ratios (a,b) $\lambda =0.45$ and (c,d) $\lambda =0.225$. Dashed boxes in (a,c) represent the corresponding locations of SEM images in (b,d), respectively. The SEM images were obtained from samples rotated $40^\circ$ about the $\tilde {z}$ axis as seen in (f) (scale bars: (b) $50\ \mathrm {\mu }\textrm {m}$ and (d) $100\ \mathrm {\mu }\textrm {m}$). Regimes: (I) meniscus deformation, (II) meniscus recession, (III) corner flow and (IV) corner transition. (e) A SEM image of cured NOA74 in a channel with $\lambda =0.225$ depicting the pinned contact line at the top of the channel sidewall. (f) Schematic of sample orientation for SEM visualization.

2.3. Free-surface morphology

We investigate the effect of channel aspect ratio $\lambda$ on the free-surface morphology by initially examining the profilometry scans in figures 3(a) and 3(c). It is observed that for both channel aspect ratios the liquid height $\tilde {h}$ at the centre of the channel decreases down the channel length. This decrease in $\tilde {h}$ at the centre of the channel results in an increase of the free-surface curvature, causing capillary pressure gradients that drive the flow. At a certain distance down the channel, $\tilde {h}$ at the centre of the channel goes to zero and the liquid splits into two filaments along the bottom corners. The filament morphology can be seen in the SEM images in figures 3(b) and 3(d).

From figures 3(a) and 3(b), it appears that the free-surface morphology can be divided into three regimes. The first is a meniscus-deformation regime (I) (or accommodation regime) where the liquid is pinned to the top of sidewalls and the top meniscus curvature increases down the length of the channel. The second is a meniscus-recession regime (II) where the liquid depins from the top of the channel wall and the meniscus begins to recede down the channel walls. The third is a corner-flow regime (III) where the liquid splits and recedes into the corners.

Examination of figures 3(c) and 3(d) also suggests the presence of three regimes. The first is the meniscus-deformation regime (I) similar to that seen in figure 3(b). However, in this case the meniscus splits into filaments prior to the liquid depinning from the top of the channel wall (see figure 3e) so that the meniscus-recession regime (II) is absent. After the splitting of the meniscus a corner-transition regime (IV) is observed, where the liquid remains pinned to the top of the channel wall. This is followed by a corner-flow regime (III) similar to that seen in figure 3(b).

The above visualizations suggest that there is a critical channel aspect ratio $\lambda _c$ at which the free-surface morphology transitions from that seen in figures 3(a) and 3(b) to that seen in figures 3(c) and 3(d). This is in agreement with experimental observations of Seemann et al. (2005). In their study, polystyrene droplets were deposited on grooves with rectangular cross-sections via vapour condensation. The polystyrene droplets flowed in the grooves and were solidified by lowering the temperature of the polymer below its glass transition temperature. The solidified samples were then characterized using atomic force microscopy.

The expression that Seemann et al. (2005) used for $\lambda _c$ was

(2.1)\begin{equation} \lambda_c=\frac{1-\sin\theta_0}{2\cos\theta_0}, \end{equation}

by assuming a circular upper meniscus contacting the bottom of the rectangular channel while being attached to the top of the channel walls. For NOA74 with $\theta _0=10^\circ$ we obtain $\lambda _c=0.42$ from (2.1), which is consistent with the free-surface morphology transition observed in figure 3. While the corner-transition (IV) regime has been previously observed experimentally, it has not, to the best of our knowledge, been accounted for in theoretical studies.

In the following sections we evaluate the importance of the free-surface morphology in model predictions. This is achieved by comparing two theoretical models describing capillary flow, where one accounts for the complexity of the free-surface morphology (§ 3.2) whereas the other assumes that it is flat (§ 3.1).

3. Mathematical modelling

Here, we describe two mathematical models for capillary flow of a non-volatile, isothermal Newtonian liquid in an open rectangular channel in contact with an ambient passive gas. We consider a liquid of density $\rho$, viscosity $\mu$, surface tension $\sigma$ and equilibrium contact angle $\theta _0$. The open rectangular channel has width $W$, height $H$ and length $L$. In this work, we use the notation $\tilde {f}$ to denote the dimensional version of a variable $f$.

3.1. Modified Lucas–Washburn model

In this model, the flow is assumed to be driven by the capillary pressure gradient caused by the circular-arc meniscus front, while viscous forces resist the flow and inertial and gravitational forces are neglected. The capillary driving force is obtained by assuming a flat upper liquid–air interface and a circular-arc meniscus front governed by fluid statics, while the viscous force is obtained by assuming a fully developed parallel flow. Through conservation of linear momentum in the axial direction an analytical expression is obtained for the position of the meniscus front $\tilde {z}_m$ as a function of time $\tilde {t}$, which in dimensional form is

(3.1)\begin{equation} \tilde{z}_m = \sqrt{\tilde{k} \tilde{t}}, \quad \mathrm{where} \ \tilde{k}=\frac{2\sigma H \zeta_o(\lambda)}{3\mu}[\cos \theta_0 (1+2\lambda)-1]. \end{equation}

Here, $\tilde {k}$ is the mobility parameter and has units of (length)$^2$/time, $\lambda =H/W$ is the channel aspect ratio and $\zeta _o(\lambda )$ is an aspect-ratio function defined as

(3.2)\begin{equation} \zeta_o(\lambda) = \frac{24}{{\rm \pi}^4 \lambda^2} \sum^{\infty}_{n=0} \frac{1}{(2n+1)^4} \left[1-\frac{\tanh[(2n+1){\rm \pi} \lambda]}{(2n+1){\rm \pi} \lambda}\right]. \end{equation}

Detailed derivation of $\zeta _o(\lambda )$ can be found in the work of Ouali et al. (2013). Equation (3.1) will be referred to as the MLW model and has been used in several other studies (Baret et al. 2007; Yang et al. 2011; Ouali et al. 2013; Sowers et al. 2016; Kolliopoulos et al. 2019).

3.2. Lubrication-theory-based model

We develop a model describing capillary flow in open rectangular channels that accounts for the non-flat shape of the upper liquid–air interface, which results in capillary pressure gradients that drive flow. This is a more complex model than the MLW model (3.1), which assumes a flat upper liquid–air interface.

3.2.1. Model geometry

We begin by considering flow in an open rectangular channel as depicted in figure 4, motivated by the experiments in § 2.3. Recall that $\lambda _c$ is the aspect ratio at which the circular upper meniscus contacts the bottom of the rectangular channel while being attached to the top of the channel sidewalls with a contact angle of $\theta _0$.

Figure 4. Schematic of liquid undergoing capillary flow in an open rectangular channel for aspect ratios (a) $\lambda \ge \lambda _c$ and (b) $\lambda <\lambda _c$.

For $\lambda \ge \lambda _c$ (figure 4a) the free-surface morphology is divided into three regimes along the $\tilde {z}$ axis as discussed in § 2.3: a meniscus-deformation regime $[0,\tilde {z}_d(\tilde {t})]$, a meniscus-recession regime $[\tilde {z}_d(\tilde {t}),\tilde {z}_m(\tilde {t})]$ and a corner-flow regime $[\tilde {z}_m(\tilde {t}),\tilde {z}_t(\tilde {t})]$. In the meniscus-deformation regime, the liquid is pinned to the top of the channel wall ($\tilde {a}=H$). The channel inlet is assumed to be fully filled ($\theta (0,\tilde {t})={\rm \pi} /2$) and the upper meniscus curvature increases down the channel length ($\theta (\tilde {z},\tilde {t})$ decreases) until the contact angle $\theta (\tilde {z}_d,\tilde {t})=\theta _0$, where $\tilde {z}_d$ is the end of meniscus-deformation regime.

We then transition to the meniscus-recession regime, where the contact angle $\theta =\theta _0$ and the liquid height starts to recede down the channel sidewalls ($\tilde {a}(\tilde {z},\tilde {t})$ decreases) until the upper meniscus contacts the channel bottom. From (2.1), here $\tilde {a}(\tilde {z}_m,\tilde {t})=W\lambda _c$, where $\tilde {z}_m$ is the meniscus position. This results in a morphology transition where the flow splits into the channel corners, leading to the corner-flow regime. In the corner-flow regime, $\theta =\theta _0$ at the channel bottom and the sidewall, and the liquid height on the sidewall $\tilde {a}(\tilde {z},\tilde {t})$ decreases from $\tilde {a}(\tilde {z}_m,\tilde {t})=W\lambda _c$ to $\tilde {a}(\tilde {z}_t,\tilde {t})=0$, where $\tilde {z}_t$ is the finger tip position.

For $\lambda <\lambda _c$ (figure 4b) the free-surface morphology is also divided into three regimes: a meniscus-deformation regime $[0,\tilde {z}_m(\tilde {t})]$, a corner-transition regime $[\tilde {z}_m(\tilde {t}),\tilde {z}_c(\tilde {t})]$ and a corner-flow regime $[\tilde {z}_c(\tilde {t}),\tilde {z}_t(\tilde {t})]$ as discussed in § 2.3. In the meniscus-deformation regime the liquid is pinned to the top of the channel wall ($\tilde {a} =H$). The channel inlet is assumed to be fully filled ($\theta (0,\tilde {t})={\rm \pi} /2$) and the upper meniscus curvature increases down the channel length ($\theta (\tilde {z},\tilde {t})$ decreases) until the contact angle $\theta (\tilde {z}_m,\tilde {t})=\theta _C$, where $\theta _C$ is the contact angle at the channel sidewall when the upper meniscus touches the channel bottom.

After the upper meniscus contacts the channel bottom, it splits into the channel corners, leading to the corner-transition regime. In this regime the liquid remains pinned to the top of the channel wall ($\tilde {a} =H$), and we assume that the contact angle at the channel bottom reaches $\theta _0$ instantaneously. To conserve mass, the contact angle at the sidewall must change from $\theta _C$ to $\theta (\tilde {z}_m,\tilde {t})=\theta _T$, where $\theta _T$ is defined in § 3.2.4 (note that $\theta _T=\theta _C$ if $\theta _0=0$). The upper meniscus curvature increases down the channel length ($\theta (\tilde {z},\tilde {t})$ at the sidewall decreases) until $\theta (\tilde {z}_c,\tilde {t})=\theta _0$, where $\tilde {z}_c$ is the finger depinning position (the position at which the liquid depins from the top of the channel wall). Once $\theta =\theta _0$ on the channel bottom and sidewall, the morphology transitions to the corner-flow regime where the liquid height on the channel sidewall $\tilde {a}(\tilde {z},\tilde {t})$ decreases from $\tilde {a}(\tilde {z}_c,\tilde {t})=H$ to $\tilde {a}(\tilde {z}_t,\tilde {t})=0$.

In the following sections we develop a mathematical model for capillary flow considering both $\lambda \ge \lambda _c$ (figure 4a) and $\lambda <\lambda _c$ (figure 4b) and accounting for the complex upper liquid–air interface morphology.

3.2.2. Governing equations

We consider mass and momentum conservation of an incompressible Newtonian liquid with constant density, given by

(3.3a)\begin{gather} \tilde{\boldsymbol{\nabla}} \boldsymbol{\cdot} \tilde{\boldsymbol{u}}=0, \end{gather}

(3.3b)\begin{gather}\rho \left[\frac{\partial \tilde{\boldsymbol{u}}}{\partial \tilde{t}}+ (\tilde{\boldsymbol{u}} \boldsymbol{\cdot} \tilde{\boldsymbol{\nabla}})\tilde{\boldsymbol{u}}\right]=-\tilde{\boldsymbol{\nabla}}\tilde{p} +\mu \tilde{\nabla}^2 \tilde{\boldsymbol{u}} + \rho \boldsymbol{\tilde{g}}, \end{gather}

where $\tilde {\boldsymbol {u}}=(\tilde {u},\tilde {v},\tilde {w})$ is the velocity field in Cartesian coordinates, $\tilde {p}$ is the liquid pressure and ${\tilde {\boldsymbol {g}}}=(\tilde {g}_x,\tilde {g}_y,\tilde {g}_z)$ is the gravitational acceleration. The no-slip and no-penetration conditions are applied along the solid walls as

(3.4)\begin{equation} \tilde{\boldsymbol{u}}=0. \end{equation}

The boundary conditions for the normal and tangential stresses at the liquid–air interface $\tilde {h}(\tilde {x},\tilde {z},\tilde {t})$ are given by

(3.5a)\begin{gather} {[\![\boldsymbol{n}\boldsymbol{\cdot} \tilde{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]}=\sigma ( \tilde{\boldsymbol{\nabla}}_s \boldsymbol{\cdot} \boldsymbol{n}), \end{gather}

(3.5b)\begin{gather}{[\![\boldsymbol{t}_1\boldsymbol{\cdot} \tilde{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]}=0, \end{gather}

(3.5c)\begin{gather}{[\![\boldsymbol{t}_2\boldsymbol{\cdot} \tilde{\boldsymbol{T}}\boldsymbol{\cdot}\boldsymbol{n}]\!]}=0. \end{gather}

Here, $\tilde {{\boldsymbol {T}}}=-\tilde {p}{\boldsymbol {I}}+\mu [{\tilde {\boldsymbol {\nabla }}} \tilde {\boldsymbol {u}}+({\tilde {\boldsymbol {\nabla }}}\tilde {\boldsymbol {u}})^\textrm {T}]$ is the stress tensor, ${\boldsymbol {I}}$ is the identity tensor, ${\tilde {\boldsymbol {\nabla }}}_s={\tilde {\boldsymbol {\nabla }}}-{\boldsymbol {n}}({\boldsymbol {n}} \boldsymbol {\cdot }{\tilde {\boldsymbol {\nabla }}})$ is the surface gradient operator, ${\boldsymbol {n}}$ is the unit outward normal vector and ${\boldsymbol {t}}_1$ and ${\boldsymbol {t}}_2$ are the two tangent vectors at the interface in the transverse and axial directions, respectively (expressions for these vectors can be found in the supplementary material available at https://doi.org/10.1017/jfm.2020.986).

Equations (3.3a) and (3.3b) are rendered dimensionless using the following scalings:

(3.6)\begin{equation} \left.\begin{gathered} (\tilde{x},\tilde{y},\tilde{z})=(Hx,Hy,Lz), \quad \tilde{t}=\frac{L}{U}t, \quad \tilde{p}=\frac{\mu U}{\epsilon H}p,\\ (\tilde{u},\tilde{v},\tilde{w})=(\epsilon Uu,\epsilon Uv,Uw), \quad \epsilon=\frac{H}{L}, \quad U=\frac{2\epsilon\sigma}{\mu}. \end{gathered}\right\} \end{equation}

Additionally, the gravitational acceleration vector is scaled as $(\tilde {g}_x,\tilde {g}_y,\tilde {g}_z)=(gg_x,gg_y,gg_z)$, where $g$ is the magnitude of the gravitational acceleration. The dimensionless parameters that arise are the Reynolds number $Re=\rho U H/\mu$ (ratio of inertial to viscous forces), the capillary number $Ca = \mu U/\epsilon \sigma$ (ratio of viscous to surface tension forces) and the Bond number $Bo=\rho g H^2/\sigma$ (ratio of gravitational to surface tension forces).

In the limits where $\epsilon ^2\ll 1$, $\epsilon Re\ll 1$ and $Bo/Ca\ll \epsilon$, the governing equations reduce to

(3.7a)\begin{gather} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z}=0, \end{gather}

(3.7b)\begin{gather}\frac{\partial p}{\partial x} = \frac{\partial p}{\partial y}=0, \end{gather}

(3.7c)\begin{gather}\frac{\partial p}{\partial z}= \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2}. \end{gather}

The boundary conditions for the normal (3.5a), transverse tangential (3.5b) and axial tangential (3.5c) stresses at the free surface reduce to

(3.8a)\begin{gather} p=-Ca^{-1} \frac{\partial_{x}^2h}{[1+(\partial_{x}h)^2]^{3/2}}=-Ca^{-1} \left[\frac{\partial_{x}h}{[1+(\partial_{x}h)^2]^{1/2}}\right]_x, \end{gather}

(3.8b)\begin{gather}0= [1-(\partial_{x}h)^2]\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)+2\partial_{x}h \left(-\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\right)-\partial_{z}h\left(\frac{\partial w}{\partial x}+\partial_xh \frac{\partial w}{\partial y}\right), \end{gather}

(3.8c)\begin{gather}0=\frac{\partial w}{\partial y}-\partial_xh \frac{\partial w}{\partial x}. \end{gather}

The normal stress balance in (3.5a) has as a special case the Young–Laplace equation $p=-Ca^{-1}\kappa$, where $\kappa$ accounts for both transverse and axial curvature contributions. However, in the limit $\epsilon ^2\ll 1$, axial curvature contributions are negligible and only the leading-order transverse curvature contributions are accounted for in (3.8a). Based on (3.7b) the $\textit{O}(1)$ term in $p$ is only dependent on $z$ and $t$, and thus the leading-order curvature term (term in brackets on the far right of (3.8a)) is actually independent of $x$ and must only depend on $z$ and $t$. The derivation of (3.7) and (3.8) can also be seen in Yang & Homsy (2006) and White & Troian (2019), who considered V-shaped channel cross-sections.

Up to this point no assumption has been made regarding the channel cross-sectional geometry. Here, we consider two geometries for the channel cross-section: (a) rectangular (figure 5a) and (b) V-shaped (figure 5b). Using these two geometries we can describe all the liquid cross-sections in figures 4(a) and 4(b) in terms of the liquid height on the solid wall $a(z,t)$ and the contact angle $\theta (z,t)$. The meniscus-deformation ($a=1$) and meniscus-recession ($\theta =\theta _0$) regimes are described using the rectangular cross-section, while the corner-transition ($a=1$) and corner-flow ($\theta =\theta _0$) regimes are described using the V-shaped cross-section.

Figure 5. Cross-sectional schematics of (a) meniscus-deformation ($a=1$) and meniscus-recession ($\theta =\theta _0$) regimes, and (b) corner-transition ($a=1$) and corner-flow ($\theta =\theta _0$) regimes.

Each cross-sectional geometry requires three additional boundary conditions to obtain expressions for $p(z,t)$ and $h(x,z,t)$: the contact-line location on the solid wall, a symmetry condition and the definition of the contact angle $\theta$. Expressions for these boundary conditions can be found in the supplementary material. We obtain expressions for $p(z,t)$ and $h(x,z,t)$ as a function of $a(z,t)$ and $\theta (z,t)$ for each regime in figure 4 by integrating (3.8a) twice with respect to $x$ and imposing the boundary conditions. The resulting $\textit{O}(1)$ expressions are

(3.9a)\begin{gather} \left.\begin{gathered} p=-\lambda\cos\theta(z,t)\\ h=1+\dfrac{\tan\theta(z,t)}{2\lambda}- \left[\dfrac{1}{4\lambda^2\cos^2\theta(z,t)}-x^2\right]^{1/2} \end{gathered}\right\} \quad \text{meniscus deformation},\end{gather}

(3.9b)\begin{gather} \left.\begin{gathered} p=-\lambda\cos\theta_0 \\ h=a(z,t)+ \dfrac{\tan\theta_0}{2\lambda}-\left[\dfrac{1}{4\lambda^2\cos^2\theta_0}-x^2\right]^{1/2} \end{gathered}\right\} \quad \text{meniscus recession},\end{gather}

(3.9c)\begin{gather} \left.\begin{gathered} p=-\dfrac{\cos\theta_0-\sin\theta(z,t)}{2} \\ h=\dfrac{\cos\theta(z,t)}{\cos(\theta(z,t)+\beta)}-\left[ \left(\dfrac{\sin\beta}{\cos(\theta(z,t)+\beta)}\right)^2-x^2\right]^{1/2} \end{gathered}\right\} \quad \text{corner transition},\end{gather}

(3.9d)\begin{gather} \left.\begin{gathered} p=-\dfrac{\cos\theta_0-\sin\theta_0}{2a(z,t)} \\ h= \dfrac{a(z,t)\cos\theta_0}{\cos(\theta_0+{\rm \pi}/4)}-\left[\left( \dfrac{a(z,t)\sin{\rm \pi}/4}{\cos(\theta_0+{\rm \pi}/4)}\right)^2-x^2\right]^{1/2} \end{gathered}\right\} \quad \text{corner flow},\end{gather}

where $\theta _0$ is the equilibrium contact angle, $\beta =\arctan (\cos \theta /\cos \theta _0)$ (see supplementary material for further details) and $\lambda$ is the channel aspect ratio. Equations (3.9a) and (3.9b) were also used by Tchikanda et al. (2004) and Nilson et al. (2006). A similar expression to (3.9c) can be found in Weislogel & Nardin (2005). The expressions in (3.9d) were also used by Romero & Yost (1996), Weislogel & Lichter (1998), Nilson et al. (2006), Yang & Homsy (2006) and White & Troian (2019). We note that to reconstruct free-surface profiles, the height profiles $h$ in (3.9c) and (3.9d) corresponding to figure 5(b) must be rotated by angles $\beta = \arctan (\cos \theta /\cos \theta _0)$ and $\beta ={\rm \pi} /4$, respectively, to match the orientation of the channel cross-section in figures 4(a) and 4(b).

3.2.3. Diffusion equations

Lenormand & Zarcone (1984) derived the following expression from system (3.7), relating the gradient in the dimensionless flux $q$ to the time derivative of the dimensionless liquid cross-sectional area $A$:

(3.10)\begin{equation} \frac{\partial A}{\partial t}=-\frac{\partial q}{\partial z}. \end{equation}

The dimensionless flux is defined as

(3.11)\begin{equation} q=\int_{A}w \, \mathrm{d}A=-\frac{\partial p}{\partial z}\bar{w}_iA, \end{equation}

where $\bar {w}_i$ is a rescaled cross-sectional-averaged dimensionless velocity. Here, $i$ is equal to $D$, $T$ or $C$ for the meniscus-deformation, corner-transition and corner-flow regimes, respectively. (As discussed below, the meniscus-recession regime is neglected.) Details of the calculation of $\bar {w}_i$ are discussed in the supplementary material.

It is evident from (3.9a)–(3.9d) that the dimensionless streamwise flux $q$ in (3.11) is a function of either the dimensionless liquid height $a(z,t)$ on the sidewall or the liquid contact angle $\theta (z,t)$ depending on the regime. Rather than considering $a(z,t)$ and $\theta (z,t)$ separately, we introduce the liquid saturation $s=\tilde {A}(a,\theta )/HW=\lambda A$ (ratio of channel cross-sectional area filled with liquid to total channel cross-sectional area). For each regime the liquid saturation $s$ is given by

(3.12a)\begin{gather} s=\frac{1}{2\lambda}\left[2\lambda-\frac{\arcsin(\cos\theta)}{2\cos^2\theta}+\frac{1}{2}\tan\theta\right], \quad \text{meniscus deformation}, \end{gather}

(3.12b)\begin{gather}s=\frac{1}{2\lambda}\left[2\lambda a-\frac{\arcsin(\cos\theta_0)}{2\cos^2\theta_0}+\frac{1}{2}\tan\theta_0\right], \quad \text{meniscus recession}, \end{gather}

(3.12c)\begin{gather}s=\frac{\lambda\hat{B}(\theta,\theta_0)}{(\cos\theta_0-\sin\theta)^2}, \quad \text{corner transition}, \end{gather}

(3.12d)\begin{gather}s=a^2\frac{2\lambda\hat{A}(\theta_0)}{(\cos\theta_0-\sin\theta_0)^2}, \quad\text{corner flow}, \end{gather}

where the geometric functions $\hat {B}$ and $\hat {A}$ can be found, respectively, in (A1) and (A2) in appendix A. Equations (3.12a), (3.12b) and (3.12d) are equivalent to expressions reported by Nilson et al. (2006).

Since the pressure $p$ in (3.9b) is constant in the meniscus-recession regime, the flux $q=0$ for this regime based on (3.11). This is because the transverse curvature gradients are zero and the only contribution to $q$ is from the $\textit{O}(\epsilon ^2)$ axial curvature gradients, for which we did not account. Nilson et al. (2006) estimated that the meniscus-recession regime size $\delta \approx L(\epsilon ^2 \lambda /2)^{1/3}$. For the microchannel dimensions considered in our study $\delta \approx 180\text {--}320\ \mathrm {\mu }\textrm {m}$, which is negligible considering the channel length is 30 mm. This estimate for $\delta$ agrees with the observations in figure 3(b) where the meniscus-recession regime size is less than $50\ \mathrm {\mu }\textrm {m}$. Therefore, effects of the meniscus-recession regime will be neglected (i.e. $\tilde {z}_d=\tilde {z}_m$ in figure 4a) and the regime transition from meniscus deformation to corner flow (for $\lambda /\lambda _c>1$) will be treated as a saturation jump.

By using (3.12a)–(3.12d) in (3.10) we obtain the following system of nonlinear partial differential equations governing the liquid saturation:

(3.13a)\begin{gather} \frac{\partial s}{\partial t}=\frac{\partial}{\partial z}\left(D_D s\frac{\partial s}{\partial z}\right), \quad \text{meniscus deformation}, \end{gather}

(3.13b)\begin{gather}\frac{\partial s}{\partial t}=\frac{\partial}{\partial z}\left(D_T s^{1/2}\frac{\partial s}{\partial z}\right), \quad \text{corner transition}, \end{gather}

(3.13c)\begin{gather}\frac{\partial s}{\partial t}=\frac{\partial}{\partial z}\left(D_C s^{1/2}\frac{\partial s}{\partial z}\right), \quad \text{corner flow}, \end{gather}

where

(3.14a)\begin{gather} D_D=\bar{w}_D\left(\frac{2\lambda^2\sin\theta\cos^2\theta}{1-\tan\theta\arcsin(\cos\theta)}\right), \end{gather}

(3.14b)\begin{gather}D_T=\bar{w}_T\frac{(\cos\theta_0-\sin\theta)^2}{4}\left(\frac{1}{\lambda\hat{B}(\theta,\theta_0)}\right)^{1/2}\left(\frac{\hat{B}(\theta,\theta_0)}{\hat{B}(\theta,\theta_0)-\tan\theta(\cos\theta_0-\sin\theta)^2}\right), \end{gather}

(3.14c)\begin{gather}D_C=\bar{w}_C \frac{(\cos\theta_0-\sin\theta_0)^2}{4}\left(\frac{1}{2\lambda\hat{A}(\theta_0)}\right)^{1/2}. \end{gather}

The quantities $D_D$, $D_T$ and $D_C$ can be thought of as dimensionless diffusion coefficients describing the interface growth.

Recall that $\lambda _c$ (see (2.1)) is the aspect ratio at which the circular upper meniscus contacts the bottom of the rectangular channel while being attached to the top of the channel sidewalls with a contact angle of $\theta _0$. When $\lambda \ge \lambda _c$ (figure 4a), the bounds of the meniscus-deformation and corner-flow regimes are ($0,z_m$) and ($z_m,z_t$), respectively, where $z_m$ is the meniscus position and $z_t$ is the finger tip position. When $\lambda <\lambda _c$ (figure 4b), the bounds of the meniscus-deformation, corner-transition and corner-flow regimes are ($0,z_m$), ($z_m,z_c$) and ($z_c,z_t$), respectively, where $z_c$ is the finger depinning position.

3.2.4. Similarity transformation

We exploit the self-similar nature of the nonlinear diffusion equations (3.13) by introducing the variable $\eta =z/\sqrt {t}$ (Romero & Yost 1996; Weislogel & Lichter 1998; Chen, Weislogel & Nardin 2006; White & Troian 2019). For $\lambda \ge \lambda _c$ shown in figure 4(a), the self-similar governing equations are

(3.15a)\begin{gather} -\frac{1}{2}\eta\frac{\textrm{d}s}{\textrm{d}\eta}= \frac{\textrm{d}}{\textrm{d}\eta}\left(D_D s\frac{\textrm{d} s}{\textrm{d}\eta}\right), \quad \eta\in[0,\delta_0\eta_0), \quad \text{meniscus deformation}, \end{gather}

(3.15b)\begin{gather}-\frac{1}{2}\eta\frac{\textrm{d}s}{\textrm{d}\eta}= \frac{\textrm{d}}{\textrm{d}\eta}\left(D_Cs^{1/2}\frac{\textrm{d} s}{\textrm{d}\eta}\right), \quad \eta\in(\delta_0\eta_0,\eta_0], \quad \text{corner flow}, \end{gather}

subject to

(3.15c)\begin{equation} s(0)=1, \quad s(\delta_0\eta_0)^{-}=s_D, \quad s(\delta_0\eta_0)^{+}=s_C, \quad s(\eta_0)=0, \end{equation}

where $\eta _0=z_t/\sqrt {t}$ is the rescaled finger tip position and $\delta _0\eta _0=z_m/\sqrt {t}$ is the rescaled meniscus position. The channel cross-section at the inlet is assumed to be fully filled and $\theta ={\rm \pi} /2$. At the end of the meniscus-deformation regime $\theta =\theta _0$, which is used in (3.12a) to calculate $s_D$.

At the beginning of the corner-flow regime $a=\lambda _c/\lambda$, which is used in (3.12d) to determine $s_C$. (Recall from § 3.2.1 that the corner-flow regime begins when $\tilde {a} = W\lambda _c$, which in dimensionless form is $a=\lambda _c/\lambda$.) Finally at the finger tip, the liquid height goes to zero. Note that $s_D=s_C$ only for $\lambda =\lambda _c$ (meniscus contacts the channel bottom at end of the meniscus-deformation regime). Equation (2.1) is used to determine $\lambda _c$, which depends only on $\theta _0$. For all $\lambda$, it is assumed in the corner-flow regime that the contact angle on the channel sidewall and bottom is always $\theta _0$, and thus independent of speed. This is the simplest assumption and allows us to focus on the influence of other problem parameters.

Two additional conditions are required to determine $\eta _0$ and $\delta _0$, which specify the bounds of each regime. The first condition is the flux matching condition given by

(3.15d)\begin{equation} \left[-D_D s\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{-}-\left[-D_C s^{1/2}\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{+}=0, \quad \mathrm{at} \ \eta=\delta_0\eta_0, \end{equation}

where the second term in each bracket accounts for the potential discontinuity in $s$ due to transitioning from the meniscus-deformation regime to the corner-flow regime. A derivation of (3.15d) can be seen in the appendix B. The second condition is that the flux approaches zero at the finger tip (i.e. $D_Cs^{1/2}\,\textrm {d}s/\textrm {d}\eta \rightarrow 0$, as $\eta \rightarrow \eta _0$). Following Romero & Yost (1996) and using (3.15b), it can be shown that to satisfy this condition, the following must be true:

(3.15e)\begin{equation} -D_Cs^{-1/2}\frac{\textrm{d}s}{\textrm{d}\eta}=\frac{1}{2}\eta_0, \quad \mathrm{at} \ \eta=\eta_0. \end{equation}

For $\lambda < \lambda _c$ shown in figure 4(b), the self-similar governing equations are

(3.16a)\begin{gather} -\frac{1}{2}\eta\frac{\textrm{d}s}{\textrm{d}\eta}= \frac{\textrm{d}}{\textrm{d}\eta}\left(D_D s\frac{\textrm{d} s}{\textrm{d}\eta}\right), \quad \eta \in[0,\delta_0\eta_0), \quad \text{meniscus deformation}, \end{gather}

(3.16b)\begin{gather}-\frac{1}{2}\eta\frac{\textrm{d}s}{\textrm{d}\eta}=\frac{\textrm{d}}{\textrm{d}\eta}\left(D_T s^{1/2} \frac{\textrm{d} s}{\textrm{d}\eta}\right), \quad \eta\in(\delta_0\eta_0,\delta_1\eta_0), \quad \text{corner transition}, \end{gather}

(3.16c)\begin{gather}-\frac{1}{2}\eta\frac{\textrm{d}s}{\textrm{d}\eta}=\frac{\textrm{d}}{\textrm{d}\eta}\left(D_C s^{1/2} \frac{\textrm{d} s}{\textrm{d}\eta}\right), \quad \eta\in(\delta_0\eta_1,\eta_0], \quad \text{corner flow}, \end{gather}

subject to

(3.16d)\begin{align} s(0)=1, \quad s(\delta_0\eta_0)^{-}=s_D, \quad s(\delta_0\eta_0)^{+}=s_T, \quad s(\delta_1\eta_0)^{-}=s(\delta_1\eta_0)^{+}=s_C, \quad s(\eta_0)\!=\!0, \end{align}

where $\delta _1\eta _0=z_c/\sqrt {t}$ is the rescaled finger depinning position (§ 3.2.1). The channel cross-section at the inlet is assumed to be fully filled and $\theta ={\rm \pi} /2$. At the end of the meniscus-deformation regime $\theta =\theta _C$ (critical angle at which the upper meniscus touches the channel bottom, calculated from $\lambda =(1-\sin \theta _C)/2\cos \theta _C$), which is used in (3.12a) to calculate $s_D$, and the contact angle at the channel bottom is $\theta =0$.

At the transition from the meniscus-recession regime to the corner-transition regime, the liquid remains pinned to the top of the channel sidewall and the upper meniscus contacts the channel bottom with the flow splitting into the channel corners. At the beginning of the corner-transition regime we assume the liquid instantaneously attains $\theta _0$ at the channel bottom and $\theta _C \rightarrow \theta _T$ at the channel sidewall. To conserve mass, we equate the amount of liquid in the channel cross-section on each side of this transition. This specifies $\theta _T$, which is calculated (via Newton's method) by setting (3.12c) equal to $s_D$. If the calculated $\theta _T\le {\rm \pi}/4$, then $s_T=s_D$. If $\theta _T>{\rm \pi} /4$, then $D_T<0$ which makes the problem ill-posed (Romero & Yost 1996). In this case, we set $\theta _T={\rm \pi} /4$, leading to a saturation jump. Equation (3.12c) is then used to determine $s_T$ based on $\theta _T$. Note that in the corner-transition regime the contact line at the channel sidewall is assumed to be pinned while the contact line at the channel bottom is allowed to move with constant contact angle $\theta _0$. At the end of the corner-transition regime and the beginning of the corner-flow regime $a=1$, which is used in (3.12d) to determine $s_C$. (Recall from § 3.2.1 that the corner-flow regime begins when $\tilde {a}=H$, which in dimensionless form is $a=1$.) Finally at the finger tip, the liquid height goes to zero.

Three additional conditions are required to determine $\eta _0$, $\delta _0$ and $\delta _1$, which are

(3.16e)\begin{gather} \left[-D_D s\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{-}-\left[-D_T s^{1/2}\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{+}=0, \quad \mathrm{at} \ \eta=\delta_0\eta_0, \end{gather}

(3.16f)\begin{gather}\left[-D_T s^{1/2}\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{-}- \left[-D_C s^{1/2}\frac{\textrm{d} s}{\textrm{d}\eta}-\frac{1}{2}s\eta\right]_{+}=0, \quad \mathrm{at} \ \eta=\delta_1\eta_0, \end{gather}

(3.16g)\begin{gather}-D_Cs^{-1/2}\frac{\textrm{d}s}{\textrm{d}\eta}=\frac{1}{2}\eta_0, \quad \mathrm{at} \ \eta=\eta_0, \end{gather}

where (3.16a) and (3.16b) are flux matching conditions (see appendix B) and (3.16c) is the condition setting the flux to zero at the finger tip.

For $\lambda \ge \lambda _c$ the system of governing equations is (3.15), whereas for $\lambda <\lambda _c$ the system consists of (3.16). What is required to solve these systems are the cross-sectional-averaged dimensionless velocities $\bar {w}_D(s)$, $\bar {w}_T(s)$ and $\bar {w}_C(s)$, which influence the values of $D_D$, $D_T$ and $D_C$ through (3.14). The cross-sectional-averaged dimensionless velocities are calculated for a given cross-section by solving (3.7c) subject to no-slip and no-penetration conditions along the solid walls and no-stress condition (3.8c) at the liquid–air interface (see supplementary material for further details).

3.2.5. Numerical methods

Velocity fields (see supplementary material) are numerically solved for with a Galerkin finite-element method using quadratic basis functions. To validate our computations, our results for $\bar {w}_D(s)$, $\bar {w}_T(s)$ and $\bar {w}_C$ are compared to results from prior studies. Results for $\bar {w}_D(s)$ and $\bar {w}_T(s)$ are in agreement with results of Tchikanda et al. (2004) and Weislogel & Nardin (2005), respectively. Results for $\bar {w}_C$ agree with results of Ayyaswamy, Catton & Edwards (1974), Ransohoff & Radke (1988) and Yang & Homsy (2006). Note that these prior studies do not consider capillary flow in open rectangular channels over the range of contact angles $\theta _0$ and aspect ratios $\lambda$ examined in the present work.

Results for $\bar {w}_D(s)$ and $\bar {w}_T(s)$ from the finite-element method simulations are fitted using Chebyshev polynomials of the first kind using the least-squares method. These Chebyshev polynomials are then used in the system of equations (3.15) ($\lambda \ge \lambda _c$) and (3.16) ($\lambda <\lambda _c$) to evaluate $D_D$ and $D_T$. Since $\bar {w}_C$ does not depend on $s$, an exact expression for $D_C$ can be obtained via (3.14c). Both nonlinear systems of equations (3.15) and (3.16) are discretized using a second-order centred finite-difference method and solved using the fsolve solver in MATLAB.

4. Results and discussion

Similarity solutions for the liquid saturation profiles $s(\eta )$ and their dependence on the channel aspect ratio $\lambda$ and equilibrium contact angle $\theta _0$ are presented first (§ 4.1). Using these similarity solutions, three-dimensional liquid height profiles are obtained to highlight the complex free-surface morphology similar to that seen in figure 3. Model predictions for the evolution of the meniscus position $\tilde {z}_m(\tilde {t})$ from the lubrication-theory-based and MLW models are then compared to experimental observations (§ 4.2). Finally, lubrication-theory-based model predictions of the finger length evolution $l_f(t)=z_t(t)-z_m(t)$ are compared to experimental results (§ 4.3).

4.1. Saturation profiles

Computed similarity solutions of the liquid saturation $s(\eta )$ for different aspect ratios $\lambda$ and an equilibrium contact angle $\theta _0=10^\circ$ are shown in figure 6(a). Solutions for $\lambda >\lambda _c$ and $\lambda <\lambda _c$ are obtained solving the system of equations (3.15) and (3.16), respectively. These similarity solutions are valid for intermediate times, when channel entrance and end effects can be neglected.

Figure 6. Effect of (a) aspect ratio $\lambda$ and (b) equilibrium contact angle $\theta _0$ on liquid saturation profiles $s(\eta )$. The two filled symbols on each curve correspond to the end of the meniscus-deformation regime and the beginning of the corner-flow regime, respectively (see figures 4a and 4b). The dashed lines indicate the saturation jump due to neglecting the meniscus-recession regime when $\lambda >\lambda _c$ (see § 3.2.3). Note that $\lambda _c=0.42$ for $\theta _0=10^\circ$.

In figure 6(a) when $\lambda >\lambda _c$ (here $\lambda _c=0.42$), $s(0)=1$ corresponds to a fully filled channel cross-section. Moving down the length of the channel, $s$ decreases monotonically and at the meniscus position (filled symbol) the flow transitions from the meniscus-deformation regime to the corner-flow regime (see figure 4a). A jump in $s$ (dashed lines) is observed because we neglected the meniscus-recession regime as discussed in § 3.2.3. In the corner-flow regime $s$ continues to decrease until $s(\eta _0)=0$ at the finger tip.

From figure 6(a), the $s(\eta )$ profiles have a non-monotonic dependence on $\lambda$, suggesting that there is an optimal $\lambda$ for capillary flow. The effect of the equilibrium contact angle $\theta _0$ on $s(\eta )$ for an aspect ratio $\lambda =0.75$ is shown in figure 6(b). Decreasing $\theta _0$ results in more capillary filling. Although in figure 6(b) we consider $\lambda >\lambda _c$, the same trend is observed for $\lambda <\lambda _c$.

The non-monotonic effect of $\lambda$ on capillary flow becomes more clear in figure 7. Here, the rescaled finger tip position $\eta _0$, meniscus position $\delta _0 \eta _0$ and finger depinning position $\delta _1 \eta _0$ (defined in § 3.2.4) are presented as a function of $\lambda$ for different $\theta _0$. The shaded areas between the curves represent the sizes of the meniscus-deformation (I), corner-flow (III) and corner-transition (IV) regimes (seen in figures 4a and 4b), which depend on $\lambda$ and $\theta _0$.

Figure 7. Lubrication-theory-based model predictions of the rescaled finger tip position $\eta _0$, meniscus position $\delta _0 \eta _0$ and finger depinning position $\delta _1 \eta _0$ as a function of channel aspect ratio $\lambda$ for (a) $\theta _0=10^\circ$, (b) $\theta _0=18^\circ$, (c) $\theta _0=32^\circ$ and (d) $\theta _0=42^\circ$. The shaded areas between the curves represent the sizes of the meniscus-deformation (I), corner-flow (III) and corner-transition (IV) regimes.

We first consider results in figure 7(a), where $\theta _0=10^\circ$. When $\lambda \gg \lambda _c$ the flow is dominated by the meniscus-deformation (I) regime. With decreasing $\lambda$, the size of the corner-flow (III) regime monotonically increases. However, the size of the meniscus-deformation (I) regime increases and then decreases, with decreasing $\lambda$. When $\lambda$ drops below $\lambda _c$, the corner-transition (IV) regime appears. As $\lambda$ is further decreased, the sizes of the corner-flow (III) and corner-transition (IV) regimes increase, while the size of the meniscus-deformation (I) regime decreases. These trends are observed for the other values of $\theta _0$ considered in figure 7(b–d).

Similarity solutions for $s(\eta )$ are used to construct three-dimensional free-surface profiles. These solutions for $s(\eta )$ are used in (3.12) to determine $\theta (\eta )$ (via Newton's method) and $a(\eta )$ for each regime. The three-dimensional free-surface profiles $h$ are determined using (3.9). Since the $h$ expressions in (3.9c) and (3.9d) for the corner-transition and corner-flow regimes correspond to figure 5(b), they require rotation by angles $\beta =\arctan (\cos \theta /\cos \theta _0)$ and $\beta = {\rm \pi}/4$, respectively, to match the channel orientation seen in figure 4. By solving the system of ordinary differential equations (3.15) and (3.16) we can construct three-dimensional free-surface profiles for any time $\tilde {t}$. Dimensional free-surface profiles $\tilde {h}$ for $\lambda >\lambda _c$ and $\lambda <\lambda _c$ after 100 s of flow of NOA74 are depicted in figure 8.

Figure 8. Simulated liquid height profiles of NOA74 after $\tilde {t} = 100\ \textrm {s}$ of flow for (a) $\lambda = 0.45$ $(\lambda >\lambda _c)$ and (b) $\lambda = 0.225$ $(\lambda <\lambda _c)$. For NOA74, $\lambda _c=0.42$.

Qualitative agreement is observed between the three-dimensional free-surface profiles in figure 8 and the profilometry measurements in figures 3(a) and 3(c). The channel is fully filled at the channel inlet and the upper meniscus bows down as we move down the length of the channel until it contacts the channel bottom and splits into the channel corners. A quantitative comparison between theory and experiment is made in the following section.

4.2. Comparison with experiments

We compare predictions of the meniscus-position time evolution $\tilde {z}_m(\tilde {t})$ from the lubrication-theory-based and MLW models to experimental observations. This comparison is made in figure 9 for the test liquids detailed in table 1. Solid and dashed lines represent the lubrication-theory-based and MLW model predictions, respectively, while experimental observations are shown as symbols. Each panel in figure 9 includes experiments for $\lambda >\lambda _c$ and $\lambda <\lambda _c$, except for the case of propylene glycol (figure 9d) where only experiments for $\lambda >\lambda _c$ were conducted.

Figure 9. Dimensional meniscus position $\tilde {z}_m$ as a function of time $\tilde {t}$ for different channel aspect ratios $\lambda$ with (a) NOA74 ($\theta _0=10^\circ$, $\lambda _c=0.42$), (b) silicone oil ($\theta _0=18^\circ$, $\lambda _c=0.36$), (c) mineral oil ($\theta _0=32^\circ$, $\lambda _c=0.28$) and (d) propylene glycol ($\theta _0=42^\circ$, $\lambda _c=0.22$). Symbols and shaded areas represent average and range of experimental results, respectively. Solid and dashed lines represent lubrication-theory-based and MLW model predictions, respectively.

When $\lambda >\lambda _c$, the lubrication-theory-based and MLW models are in good agreement with experiments for NOA74 (figure 9a) and silicone oil (figure 9b). For mineral oil (figure 9c) and propylene glycol (figure 9d), the lubrication-theory-based model agrees well with the experiments but the MLW model underpredicts $\tilde {z}_m(\tilde {t})$.

For NOA74 (figure 9a) and silicone oil (figure 9b), the lubrication-theory-based model agrees with experiments for $\lambda <\lambda _c$, but the MLW model overpredicts $\tilde {z}_m(\tilde {t})$. For mineral oil (figure 9c), both models overpredict $\tilde {z}_m(\tilde {t})$ for $\lambda <\lambda _c$. For propylene glycol (figure 9d), at larger $\lambda$ the lubrication-theory-based model predictions agree better with experiments than predictions of the MLW model. However, as $\lambda \rightarrow \lambda _c$ the MLW model predictions agree with experiments, whereas the lubrication-theory-based model overpredicts $\tilde {z}_m(\tilde {t})$.

To further compare theory and experiment we consider the dimensionless mobility parameter $k=z_m^2/t$, which can be thought of as a diffusion coefficient driving the growth of the meniscus position. In the lubrication-theory-based model,

(4.1a)\begin{equation} k = (\delta_0\eta_0)^2, \end{equation}

since $\delta _0\eta _0=z_m/\sqrt {t}$ from (3.15a) and (3.16a). Computed values of $\delta _0 \eta _0$ for different $\lambda$ and $\theta _0$ are shown in figure 7. In the MLW model,

(4.1b)\begin{equation} k=\tilde{k}/UL, \end{equation}

where $\tilde {k}$ is defined in (3.1), $U$ is the characteristic velocity (see § 3.2) and $L$ is the channel length. Experimentally, $k$ is determined by fitting the function $z_m=(kt)^n$ to experiments similar to those shown in figure 9 using nonlinear regression. A comparison of theoretically predicted and experimentally determined $k$ values as a function of $\lambda$ and $\theta _0$ is shown in figure 10.

Figure 10. Dimensionless mobility parameter $k$ as a function of aspect ratio $\lambda$ for (a) NOA74 ($\theta _0=10^\circ$, $\lambda _c=0.42$), (b) silicone oil ($\theta _0=18^\circ$, $\lambda _c=0.36$), (c) mineral oil ($\theta _0=32^\circ$, $\lambda _c=0.28$) and (d) propylene glycol ($\theta _0=42^\circ$, $\lambda _c=0.22$). Filled symbols represent experimental results, dashed and solid lines represent MLW and lubrication-theory-based model predictions, respectively, and shaded areas represent $\lambda <\lambda _c$. Insets: fitted scaling exponent $n$ (from $z_m = (kt)^n$) as a function of aspect ratio $\lambda$ for liquids with different equilibrium contact angles $\theta _0$. Solid lines represent lubrication-theory-based and MLW model predictions of $n=0.5$.

When $\lambda \gg \lambda _c$, the lubrication-theory-based and MLW model predictions are indistinguishable regardless of $\theta _0$. In this limit the finger contribution to the flow becomes negligible and the flow asymptotically approaches that between two parallel plates. Note that when $\lambda = H/W\gg \lambda _c$, the effects from the channel bottom become negligible compared with the effects from the sidewalls. As $\lambda \rightarrow 1$ the two model predictions begin to deviate because the viscous resistance in the fingers becomes significant, and fingers are not accounted for in the MLW model. However, both model predictions for $k$ are in reasonable agreement with experiments. As $\lambda \rightarrow \lambda _c$ the lubrication-theory-based model agrees better with experiments, compared to the MLW model (except in figure 10d). The disagreement between the MLW model and experiments becomes more evident as $\lambda$ decreases because the free-surface height profiles (figures 3b and 3d) become more non-uniform. Therefore, the assumption of a flat free-surface height profile used in the MLW model becomes less valid with decreasing $\lambda$.

However, the lubrication-theory-based model does not always perform better than the MLW model (figure 10d). As $\theta _0\rightarrow {\rm \pi}/4$ the transverse free-surface curvature in the corner-transition and corner-flow regimes vanishes and the flow in these regimes is then driven by axial curvature gradients (Yang & Homsy 2006). These axial curvature gradients are neglected in the lubrication-theory-based model (see (3.8a)) because they are $\textit{O}(\epsilon ^2)$ terms. When the axial curvature contributions are included in the pressure, then

(4.2)\begin{equation} p = -Ca^{-1}\left[\mathcal{C}_T(s)+\epsilon^2\mathcal{C}_A\left(s,\partial_z s,\partial_z^2 s\right)\right], \end{equation}

where $\mathcal {C}_T(s)$ and $\mathcal {C}_A(s,\partial _z s,\partial _z^2 s)$ are the transverse and axial curvature contributions, respectively. Including the axial curvature gradients in the model results in fourth-order spatial derivatives of $s$, for which we were unable to find similarity solutions. An additional reason for the deviation between the lubrication-theory-based model predictions and the experimental observations as $\theta _0\rightarrow {\rm \pi}/4$ is that the lubrication approximation itself is expected to become less accurate as the contact angle increases.

For predicting $z_m(t)$, the choice of using the lubrication-theory-based or MLW model depends on $\lambda$ and $\theta _0$. For large $\lambda$ both models give similar results, so the MLW model is preferred because of its simplicity, but for small $\lambda$ the lubrication-theory-based model is more accurate. Additionally, the lubrication-theory-based model performs better for smaller $\theta _0$ compared to the MLW model. However, it neglects key physical contributions as $\theta _0\rightarrow {\rm \pi}/4$ which results in poorer agreement. Nevertheless, when additional information is desired, such as the free-surface morphology (figure 8) or finger dynamics (§ 4.3), then the lubrication-theory-based model must be chosen.

4.3. Finger dynamics

A key advantage of the lubrication-theory-based model compared with the MLW model is that it describes finger dynamics. The dimensionless finger length $l_f(t)$ is defined as the distance between the dimensionless finger tip position $z_t(t)$ and the dimensionless meniscus position $z_m(t)$ which can be seen in figure 4. Predictions of the dimensionless finger length time evolution $l_f(t)$ by the lubrication-theory-based model are compared to experimental observations in figure 11. Model predictions for $l_f$ were computed using

(4.3)\begin{equation} l_f=z_t(t)-z_m(t)=(1-\delta_0)\eta_0 t^{1/2}, \end{equation}

where $\eta _0=z_t/\sqrt {t}$ and $\delta _0\eta _0=z_m/\sqrt {t}$; values for $\eta _0$ and $\delta _0\eta _0$ are shown in figure 7.

Figure 11. Effect of (a) aspect ratio $\lambda$ and (b) equilibrium contact angle $\theta _0$ on dimensionless finger length $l_f$ as a function of dimensionless time $t$. Filled symbols and shaded areas represent average and range of experimental results, respectively. Solid lines represent lubrication-theory-based model predictions. The liquids are mineral oil ($\theta _0 =32^{\circ }$) and silicone oil ($\theta _0 =18^{\circ }$).

The effects of $\lambda$ and $\theta _0$ on $l_f(t)$ are shown in figures 11(a) and 11(b). Solid lines and symbols represent lubrication-theory-based model predictions and experimental observations, respectively. It is important to note the shorter dimensionless times compared to the meniscus-tracking experiments (figure 9) were caused by the finger tip moving out of the field of view. In figure 11(a), experimental observations are compared to model predictions for mineral oil ($\lambda _c=0.28$), where increasing channel aspect ratio $\lambda$ results in a decrease in $l_f(t)$. Good agreement between theory and experiments is observed, with model predictions being within the range of experimental observations for all but one trial (i.e. $\lambda =0.3$ in figure 11a). A possible reason for this discrepancy is the error in identifying the position of the finger tip where the liquid height goes to zero, which is challenging during the meniscus-tracking experiments where the stage is moving (figure 2a).

In figure 11(b), experimental observations for mineral oil ($\theta _0=32^\circ$, $\lambda _c=0.28$) and silicone oil ($\theta _0=18^\circ$, $\lambda _c=0.36$) are compared to model predictions for $\lambda =0.23$ ($\lambda <\lambda _c$). Increasing $\theta _0$ leads to a decrease in $l_f(t)$. Model predictions are in good agreement with experimental observations in this case as well. Hence, the lubrication-theory-based model accurately captures effects of $\lambda$ and $\theta _0$ on $l_f(t)$. Additionally, the good agreement between theory and experiment in figure 11 emphasizes the importance of accounting for the corner-transition regime (figure 4b) when describing the finger dynamics in open rectangular channels when $\lambda <\lambda _c$.

We make two additional comments about the fingers. First, although curvature of the channel corners may affect finger length (Chen et al. 2006; Gerlach et al. 2020), the good agreement observed here suggests that the corner curvature does not play a significant role in our experiments. Second, although the lubrication-theory-based model overpredicts the meniscus position for mineral oil (figure 9d), it accurately predicts the finger length (figure 11) which is the difference between the positions of the meniscus and finger tip (4.3). This indicates that while the lubrication-theory-based model may not predict absolute positions accurately, it can predict relative positions accurately. A deeper understanding of the reasons for this likely requires a comparison of the lubrication-theory-based model and direct numerical simulations, which is beyond the scope of the present work.

5. Conclusions

In this work we combine theory and experiment to examine capillary-flow dynamics in open rectangular microchannels. For open microchannels, the free surface greatly influences the capillary-flow dynamics. We visualize the free-surface morphology and its dependence on the channel aspect ratio $\lambda$ using SEM and profilometry. The SEM images suggest a qualitative difference in the free-surface morphology at $\lambda =\lambda _c$ and highlight the significance of a corner-transition regime when $\lambda <\lambda _c$, which was not accounted for in previous studies of capillary flow in open rectangular channels.

Effects of the free-surface morphology on capillary-flow dynamics were examined using two theoretical models. The first model is a MLW model, which assumes a flat free surface and has been extensively used in prior studies. The second model is a self-similar lubrication-theory-based model, which was developed to account for the complex free-surface morphology observed in the experiments.

Predictions of the lubrication-theory-based and MLW models were compared to complementary flow visualization experiments over a wide range of channel aspect ratios $\lambda$ and equilibrium contact angles $\theta _0$. For large $\lambda$, predictions from the two models are indistinguishable, since free-surface morphology effects are negligible. However, for smaller $\lambda$ the lubrication-theory-based model is in better agreement with experiments since the influence of fingers is significant. Additionally, the lubrication-theory-based model agrees better with experiments for smaller $\theta _0$, although as $\theta _0\rightarrow {\rm \pi}/4$ the agreement worsens because important axial curvature contributions are neglected. Finally, the lubrication-theory-based model predictions accurately capture the finger length time evolution seen in experiments over a range of $\lambda$ and $\theta _0$, further highlighting the importance of accounting for the corner-transition regime in the model formulation when $\lambda <\lambda _c$.

Our lubrication-theory-based model and flow visualization experiments reveal the limitations of the widely used MLW model. In addition, our results significantly advance physical understanding of capillary-flow dynamics in open rectangular microchannels, which play a key role in a number of technological applications (§ 1) but are less well studied than V-shaped channels. Finally, the systematic development of our model allows for it to be readily extended to incorporate other important phenomena such as gravitational effects, solute transport and evaporation.