Capillary imbibition of monodisperse emulsions in confined microfluidic channels
Masoud Norouzi Darabad, Sagnik Singha, Jerzy Blawzdziewicz, Siva A. Vanapalli, Mark W. Vaughn
1. Introduction
Capillary imbibition is a process in which a fluid driven by surface forces penetrates into a small pore. This process has various applications in science and engineering including lab-on-a-chip devices
1, 2 , oil recovery , and ink-jet printing . In addition to its applications, capillary imbibition provides an unconventional method to investigate the microscale behaviour of dispersed-phase materials. By confining the dispersed-phase flow in a thin-rectangular capillary to achieve a two-dimensional-like geometry of a single droplet layer, the differences in the motion between phases is enhanced. In order to describe the behaviour of an imbibed emulsion, we must first understand the imbibition of a simple liquid. The oldest report of this phenomena is found in the manuscripts of Leonardo da Vinci in 15 th century , but the detailed dynamics were not understood until the 20 th century. In 1921, EW Washburn analysed the capillary imbibition of Newtonian fluids in a viscous dominated regime. He reported that the square of the penetration length is a linear function of imbibition time. In the early moments of imbibition, however, inertial forces may dominate, causing the dynamics to deviate from Washburn's predictions
7, 8 . Furthermore, interpretation can be complicated by flow-induced changes in the contact angle . The imbibition process is unsteady and temporal changes in the capillary number also affects, the dynamic contact angle and results in the deviations from Washburn’s model
10, 11 . For complex fluids, capillary imbibition has numerous biological and industrial applications. Examples include blood flow in microcapillaries , flow of polymers in nano-porous materials , fabrication of nano-composite materials and optofluidics . There are numerous studies for imbibition of single phase non-Newtonian fluids
14, 18-23 , however there are only a few studies of capillary imbibition of suspensions. Zhou et. al. studied the imbibition of blood suspensions and reported a failure of penetration due to capillary clogging as a result of segregation of red blood cells. Holloway et al. studied the imbibition of a concentrated suspension of rigid spherical particles and reported that the dynamics of the imbibition obey Washburn's model, however the slope term is smaller than that predicted by using the effective viscosity of the suspension. They attributed this phenomenon to the accumulation of particles behind the meniscus and to the resulting non-uniformity of concentration in the axial direction. However, their analysis applies only to systems where the particle diameter d is much smaller than the capillary diameter. There does not appear to be a systematic investigation of the imbibition of suspensions under strong confinement conditions d/h = O(1) (where h is the smallest dimension of the channel). In contrast, the pressure driven flow of suspensions in narrow channels is well studied, especially in the Hele-Shaw geometry . Pressure driven droplets/particles in Hele-Shaw cells exhibit fascinating collective behaviour such as vibrational modes in 1D ordered droplet arrays , dislocations
29, 31 fingering instabilities in 2D ordered arrays, and Burger shock-waves developing in disordered droplet monolayers . Capillary imbibition of monodisperse emulsions in confined microfluidic channels
Masoud Norouzi Darabad, a Sagnik Singha, b Jerzy Blawzdziewicz, b,c
Siva A. Vanapalli a and Mark W. Vaughn a We investigate imbibition of a monodisperse emulsion into a low-aspect ratio microfluidic channel with the height h comparable to the droplet diameter d . For confinement ratio d/h = 1.2, the tightly confined disk-like droplets in the channel move more slowly compared to the average suspension velocity. Behind the meniscus that drives the imbibition, there is a droplet-free region, separated from the suspension region by a sharp concentration front. The suspension exhibits strong droplet density and velocity fluctuations, but on average, the suspension domain remains uniform. For weaker confinement, d/h = 0.65, the spherical droplets move faster than the average suspension flow, resulting in the formation of a dynamically unstable high-concentration region near the meniscus. We describe the macroscopic suspension dynamics using linear transport equations for the particle-phase flux and suspension flux that are driven by the local pressure gradient. A dipolar particle interaction model explains the observed large density and velocity fluctuations in terms of the dynamics of elongated particle clusters with different orientations. a. Dept. of Chemical Engineering, Texas Tech University, Box 43121 Lubbock, TX 79409, USA. b. Dept. of Mechanical Engineering, Texas Tech University, Box 41021 Lubbock, TX 79409, USA. c. Dept. of Physics and astronomy, Texas Tech University, Box 41501 Lubbock, TX 79409, USA. he observed complex collective dynamics stems, primarily, from far-field dipolar Hele-Shaw hydrodynamic interactions between strongly confined particles
29, 30, 32 . The effects of higher-order contributions on the microscale evolution
29, 35 and the macroscopic dynamics have also been investigated. The presence of particles can either create or supress instabilities near the suspension front . For example, the presence of particles in a suspension flow in a radial Hele-Shaw cell can result in viscous fingering. Additionally, under strong confinements, particles can create a dense particle-band at the liquid-air interface, which can break-up due to the deformations of the interface. Here we study capillary imbibition of emulsions in thin, low aspect-ratio rectangular microfluidic channels. The droplet diameter d is comparable to the channel height h , and the droplets form a monolayer in this Hele-Shaw geometry. We consider two characteristic cases of the imbibition process, depending on the ratio between droplet diameter d and the channel height h . In the first case, d/h >1, the droplets are flattened between the channel walls, and therefore have a circular disk-like shape. In the second case, d/h <1, the droplets remain spherical during the imbibition process. We explain how microscale droplet interactions result in distinct dynamics of the imbibition at macroscale in these two cases. For disk-like droplets, we observe formation of the droplet-free zone behind the meniscus with a sharp boundary between the droplet-free region and suspension region. For spherical droplets, as the particles accumulate behind the meniscus the densely packed accumulation region undergoes fingering instability that affects particle distribution and may cause strong deformation of the meniscus. To shed light on the phenomenon observed in our capillary imbibition experiments, we use the linear transport equations to describe suspension flow in the channel, carry out microscopic analysis of particle cluster motion, and perform numerical simulations of particle motion using a dipolar-interaction model. We compare our experimental measurements with a benchmark model based on pressure driven-flow of hard-sphere monolayers in a parallel wall channel. Materials and methods
Microfluidic device fabrication
Oil-in-water emulsions were prepared using a flow-focusing device
38, 39 with orifice width and length of 20 µm and 70 µm, respectively. The silicon wafer mould containing the designed flow-focusing devices was prepared using soft lithography techniques and the resulting channels were 50 ± Emulsion preparation
Emulsions were prepared using water containing 2 wt% sodium dodecyl sulfate (SDS) (Sigma Aldrich) as the continuous phase and bromohexadecane (Sigma Aldrich) as the dispersed phase. The 2 wt% SDS corresponds to 8.6 times the critical micelle concentration . The addition of the SDS is necessary to prevent droplet coalescence. The density of the dispersed phase, 0.999 g/mL at 25 C is close to water, so droplets are neutrally buoyant. Syringe pumps (PHD 2000, Harvard Apparatus) along with 1 ml plastic syringes are used to pump the fluids to the microfluidic device. To make monodisperse oil in water emulsions, bromohexadecane is pumped through the center channel of the flow focusing device and then sheared in the orifice by SDS solution flowing from the side channels. The interfacial tension between bromohexadecane and SDS was measured to be 0.11 mN/m using the pendant drop method. Stock emulsions were prepared with two different droplet diameters: d=60 µm and d=65 µm, both with ϕ = 0.50 oil volume fraction. The emulsions were then stored in microcentrifuge tubes (Fisher Scientific) and diluted to ϕ = 0.4, 0.3, 0.2, and 0.1 by adding 2 wt% SDS solution. Experimental protocol
The system geometry is illustrated in Fig. 1a. Imbibition experiments were done using rectangular borosilicate glass capillaries (Vitrocom) that are l = 5 cm long with a width to height ratio w/h =
10. Two different sizes of capillaries are used, one with width w = 500 µm and height h = 50 µm and the other one with w = 1000 µm and h = 100 µm. Prior to conducting an experiment, the emulsions were carefully mixed by pipetting to remove non-homogeneity in droplets concentration due to slight creaming of the droplets. Then, using the pipette, an emulsion (approximately two times the volume of the capillary tube) is carefully placed in contact with the glass capillary inlet. The imbibition process starts immediately and is monitored using an inverted microscope (Nikon Eclipse- TiU) and a high-speed imaging camera (Phantom v710 12-bit, Vision Research). Videos were recorded at a maximum of 700 fps at 2X magnification with a pixel resolution of 10 µm per pixel. Experiments were recorded either at “Location 1” that covers the region from the capillary inlet to the position 1 cm away from the inlet (0 < x/l < 0.2) or "Location 2" that covers the region from 1.5 cm 2.5 cm away from the inlet (0.3 < x/l <0.5). These two locations are illustrated in the schematic of Fig. 1a.
Image processing
Image processing is implemented using a custom MATLAB (R2020a) program to determine the positions of the droplets and the meniscus. Prior to the imbibition process, a background image of the field of view is recorded. Since the droplet boundaries and the menisci appear as dark pixels, experimental frames and background image are inverted prior to background subtraction. Subtracted images are converted to a binary image with a threshold chosen to remove noise. The coordinates of the centre of the meniscus are used as its position. The results f the MATLAB program were verified by using ImageJ to compare random samples to manual measurements. To find the centroids of the droplets, the images are thresholded, and the connected regions identified using the bwconncomp function from the MATLAB (R2020a) image processing toolbox. A connected region with size comparable to a droplet and high circularity was identified as a potential droplet. The centroids of individual droplets can be used to extract their trajectories. We have created a custom MATLAB (R2020a) code that uses a Kalman filter and the Hungarian algorithm for multiple object tracking. The extracted tracks are inspected, and faulty tracks removed. Area fraction measurements
It is essential to monitor the area fraction 𝜙 𝑎 of the emulsion inside the capillary, since it can differ from the volume fraction 𝜙 of the original emulsion. These area fraction differences can stem from phenomena such as crowding of the droplets at the inlet of the capillary, the different flow rates of the droplets and the continuous phase fluid, and non-homogeneous droplet distribution in the original emulsion due to the creaming. The area fraction is measured directly by counting the number of droplets. See Fig. S1 and S2, ESI Section 1 for the evolution of the area fraction and the average area fraction for spherical and disk-like droplets. Using these data, we report an average area fraction 𝜙̅ 𝑎 for each experiment. Dimensionless numbers
The Reynolds number Re and the capillary number Ca for this system are defined as Re = ρuh/µ and Ca = µu/σ where ρ is density, µ is viscosity, σ is surface tension, and u is the superficial velocity of the continuous phase. For a capillary with height h = 50 µm, we have Re ≈
1 at location 1 and Re ≈ h = 100 µm are Re ≈ Re ≈ Ca is less than 0.01 in all experiments. According to the above values, inertial effects are relatively small, and Stokes-flow conditions approximately apply. The droplets are not deformed by the flow (except at the channel inlet for d/h <2) and are either circular or spherical in our experiments. Simulations
To elucidate the observed droplet dynamics, we performed numerical simulations of a simplified system in which the droplets interact via the leading-order Hele – Shaw dipolar scatter flow fields
29, 30, 32 , and the near-field effects are modelled using a short-range pairwise-additive repulsive potential, which prevents overlap. Previously, similar techniques were used for systems interacting via dipolar
26, 30 , quadrupolar , and combined dipolar and quadrupolar flows produced by the particles. To replicate the effect of side walls, we used a flow reflection method
45, 46 . For macroscopically homogeneous systems, periodic boundary conditions in the flow direction were applied. For inhomogeneous systems with a sharp concentration front, the leading portion of the suspension was simulated explicitly, and the trailing portion (far behind the front) was implicitly modelled by a repulsive potential acting in the flow direction and moving with the suspension velocity. Results and discussion
System description
A schematic of the experimental cell and two typical snapshots of the imbibition process are depicted in Fig. 1. In our analysis of the experiments, we use the coordinate system where the imbibition occurs in the axial direction x , and the channel width and height are along the y and z directions. The position of the meniscus at time t is denoted by L(t).
The entrance location 1 and the mid-channel location 2 observation locations are marked by the dashed red lines (Fig. 1a). We performed two sets of experiments: one set for strongly confined disk-like droplets and the other for moderately confined spherical droplets. In the experiments with strongly confined disk-like (Fig. 1b(i)), an emulsion of diameter d = 60 µm droplets is imbibed into a capillary with height h = 50 µm (confinement ratio d/h = 1.2). After entering the capillary, the drops deform to a disk-like shape because they are squeezed between the top and Fig 1. System description. (a) Schematic of the system. The glass capillaries were 5 cm long with an aspect ratio h/w = 0.1. The imbibition experiments were recorded at two different locations, location 1: x = 0 cm to 1 cm and location 2: x = 1.5 cm to 2.5 cm. L(t) is the penetration length of the meniscus at time t. (b) (i) Top-view image of capillary imbibition of disk-like droplets. Droplets with diameter d = 60 µm imbibe into a capillary with h = 50 µm and deform to a disk-like shape. The inset schematic shows the droplets in the y-z plane. Disk-like droplets move slower than the continuous phase and a droplet-free region forms behind the interface. (ii) Top-view image of capillary Imbibition of spherical droplets. Droplets with d = 65 µm imbibe into a capillary with h = 100 µm. Spherical droplets flow faster than continuous phase and form a droplet-band behind the interface. The inset schematic shows the droplets in the y-z plane. The droplet-free region and droplet-band are shown by red rectangles. ottom walls. For moderately confined spherical droplets (Fig. 1b(ii)), an emulsion of d = 65 µm droplets are imbibed into a capillary with h = 100 µm ( d/h = 0.65). The moderately confined droplets remain spherical, but they are large enough so that only a single layer forms in the capillary. For both spherical and disk-like droplets, experimental measurements were conducted using five different volume fractions, ranging from 0 to 0.5. Two trials to monitor the imbibition process at locations 1 and 2 were performed for each volume fraction. An image of a dispersion of strongly confined disk-like droplets imbibing into the capillary channel (Fig. 1b(i)) shows that two distinct regions develop during imbibition: there is a droplet-free region behind the meniscus, followed by a suspension region. This behaviour indicates that strongly confined droplets move slower than the continuous phase fluid. In contrast, Fig. 1b(ii) shows that moderately confined spherical droplets flow faster than the continuous phase. They catch up with the meniscus and create an adjacent dense droplet-band. In the following two sections, we analyse the droplet distribution and the dynamic of the imbibition process in the disk-like and spherical droplet systems. Imbibition of disk-like droplets 3.2.1 Droplet distribution and overall imbibition dynamics
An examination of strongly confined systems depicted in Fig. 1b (i) and Fig. 2 reveals several characteristic features of the droplet distribution. First, the droplet-free and suspension regions are separated by a sharp concentration front. Second, the droplet distribution in the suspension undergoes considerable local fluctuations associated with formation of local droplet clusters. Moreover, large droplet-velocity fluctuations are observed. A similar behavior was seen for all area fractions investigated. To determine the overall droplet distribution, the visible part of the channel is divided into 13 equal length sections of 500 μm , and the area fraction is measured in each section. Fig. 2 presents the result for an experiment with an average area fraction of 𝜙̅ 𝑎 =0.3. The first 12 sections (0 to 600 𝜇𝑚 ) are from the suspension region, and the last section (6
00 to 65 𝜇𝑚 ) belongs to the droplet-free region. The bottom panel shows the measured value of the area fraction for each of the sections marked in the top panel. The results indicate considerable local droplet density fluctuations, but no systematic area-fraction gradient was detected, consistent with qualitative observations of the particle distribution. To illustrate the overall imbibition dynamics, Fig. 3 shows the time 𝑡 fill required for the emulsion to fill the entire capillary tube, plotted versus the average area fraction of the emulsion measured at locations 1 and 2. For experiments recorded at location 1 (Fig. 3a), t fill is determined by measuring the time interval between the meniscus formation at the inlet of the channel and the termination of droplet motion. For experiments recorded at location 2 (Fig. 3b), the imbibition time is evaluated from the observed time interval between the meniscus arrival at location 2 and termination of droplet motion. To this value, we add the estimated time it takes the meniscus to reach the imaged region; the estimate was obtained using Washburn analysis (details in ESI Section 2). We find that the results determined from observations at locations 1 and 2 are consistent. As depicted in Fig. 3, 𝑡 fill is an increasing function of the area fraction. The imbibition at higher area fractions is slower because the presence of droplets results in increased viscous dissipation and therefore, an increased pressure drop opposing the Laplace pressure of the meniscus. In the next section, we characterize the imbibition dynamics using appropriate transport coefficients. In the effective medium approximation, the suspension flow in the channel can be determined using a set of constitutive transport equations for the average particle velocity 𝑈̅ and the volumetric suspension velocity 𝑢̅ : 𝑈̅ = −𝜈̅ 𝑡𝑝 ∇ || 𝑝 𝑠 ( -1) Fig. 2. Distribution of disk-like droplets in the channel
A portion of the channel is divided into 13 equally sized sections and the area fraction in each section is measured. The first 12 sections (from 0 to 6000 μm ) belong to the suspension region and the last section belongs to the droplet-free region. There is a sharp increase in the area fraction due to the sharp boundary between droplet-free and suspension region. The system exhibits considerable local drop density fluctuations, but no systematic area-fraction gradient was detected, consistent with qualitative observations of the particle distribution. The snapshot in this figure is taken from an experiment with 𝜙̅ 𝒂 = Fig. 3. Filling time of capillary tube for disk-like droplets.
The imbibition time for an emulsion of disk-like droplets to fill the capillary tube measured at (a) Location 1 and (b) Location 2 of the channel. The filling time is an increasing function of the area fraction because presence of the droplets increases the viscous dissipation resulting in an increase in pressure drop that opposes the Laplace pressure of the meniscus. ̅ = −𝜈̅ 𝑝𝑝 ∇ || 𝑝 𝑠 ( -2) These two equations, written here for the axial components, link the local velocities 𝑈̅ and 𝑢̅ to the local pressure gradient ∇ || 𝑝 𝑠 via two linear transport coefficients 𝜈̅ 𝑡𝑝 and 𝜈̅ 𝑝𝑝 . In addition, we introduce the transport equation 𝑈̅ = 𝜒 𝑐 𝑢̅ ( -3) that describes particle convection from the average flow. According to these equations, the particle-phase mobility coefficient 𝜈̅ 𝑡𝑝 , the effective permeability coefficient 𝜈̅ 𝑝𝑝 and the convective transport coefficient 𝜒 𝑐 are related: 𝜈̅ 𝑡𝑝 = 𝜒 𝑐 𝜈̅ 𝑝𝑝 ( -4) The transport coefficients 𝜈̅ 𝑡𝑝 , 𝜈̅ 𝑝𝑝 , and 𝜒 𝑐 depend on the confinement ratio d/h , area fraction 𝜙̅ 𝑎 , and flow-induced local drop distribution. Our goal is to estimate the parameters 𝜈̅ 𝑡𝑝 , 𝜈̅ 𝑡𝑝 and 𝜒 𝑐 based on the capillary imbibition data. Since the fluid is incompressible, the volumetric suspension velocity 𝑢̅ is the same as the average fluid velocity in the droplet-free region 𝑢̅ 𝑢̅ = 𝑢̅ ( -5) Thus, the suspension velocity 𝑢̅ can be determined from the measurements of the meniscus velocity. The incompressibility condition eqn (3- ) also allows the pressure gradient in the suspension region, ∇ || 𝑝 𝑠 , to be related to the pressure gradient in the clear fluid region, ∇ || 𝑝 𝑐 , 𝜈̅ 𝑝𝑝 ∇ || 𝑃 𝑠 = 𝜈̅ ∇ || 𝑝 𝑐 ( -6) where 𝜈̅ = 𝜂 −1 𝐻 is the permeability of the channel filled with the clear continuous-phase fluid. The total pressure drop in the channel is ∆P = −(𝑙 𝑐 ∇ || 𝑃 𝑐 + 𝑙 𝑠 ∇ || 𝑃 𝑠 ) ( -7) where 𝑙 𝑐 and 𝑙 𝑠 are the instantaneous lengths of the droplet-free and droplet-dispersion regions and ∆𝑃 = 2 𝜎cos𝜃/ℎ is estimated from the Young-Laplace equation. The contact angle 𝜃 is not directly observable in our experiments, but we estimate it from single-phase imbibition experiments using 2 wt% SDS solution. Using eqn. (3-6), and (3-7), we find the permeability coefficient of the suspension filled channel as 𝜈̅ 𝑝𝑝 𝜈̅ = 𝑙 𝑠 𝑙 − 𝑙 𝑐 ( -8) where 𝑙 = (𝜈̅ Δ𝑃)/𝑢̅ is the penetration length for which the clear fluid that moves with the velocity 𝑢̅ in a single-phase experiment. The velocity 𝑢̅ and the lengths of the clear fluid and suspension regions 𝑙 𝑐 and 𝑙 𝑠 are measured during the emulsion imbibition process (details in ESI Section 3) then used to evaluate the permeability coefficient 𝜈̅ 𝑝𝑝 . To determine the convective transport coefficient 𝜒 𝑐 we use eqn (3-3). The average droplet velocity is evaluated from the motion of the suspension front — since there is no discernible macroscopic gradient of the area fraction 𝜙 ̅ 𝑎 (Fig. 2), the front velocity is the same as the droplet-phase velocity 𝑈̅ . After evaluating the transport coefficients 𝜈̅ 𝑝𝑝 and 𝜒 𝑐 , the particle mobility 𝜈̅ 𝑡𝑝 is obtained using eqn (3-4). All transport coefficient values evaluated from subsequent video frames during a given experiment are averaged to report the final value. Fig. 4 shows the convective particle transport coefficient 𝜒 𝑐 vs the area fraction 𝜙̅ 𝑎 , evaluated from the measurements at location 1 (Fig. 4a) and location 2 (Fig. 4b). Location 1 (close to the entrance) tends to yield higher values of 𝜒 𝑐 than location 2, even at the lowest area fractions considered. The observed location dependence of 𝜒 𝑐 may stem from a decreased drop – wall lubrication friction associated with the higher droplet velocity at location 1. Other nonlinear effects not accounted for in our analysis (such as inertial effects and the evolution of particle distribution) may also affect 𝜒 𝑐 . The results in Fig. 4 show that the convective transport coefficient 𝜒 𝑐 increases with 𝜙 ̅ 𝑎 and approaches values close to unity for higher area-fractions. This is because at higher area-fractions, the particle phase produces stronger hydrodynamic resistance for the flow of the continuous-phase fluid in the space between the particles. Therefore, both the fluid and the particles translate with the average velocity close to the overall suspension velocity. This mechanism can be qualitatively modelled by approximating the particle phase as a porous medium through which the fluid phase flows with the permeability coefficient 𝜇 𝑠 . The particle medium moves with respect to the channel walls, with the friction coefficient 𝑓 . By balancing the wall friction force with force resulting from fluid resistance, we obtain the expression (Details in the ESI Section 5) Fig 4. Convective transport coefficient 𝝌 𝒄 of disk-like droplets. Convective transport coefficients χ c of disk-like droplets vs area fraction for (a) Location 1 and (b) Location c increases with 𝜙̅ 𝑎 . This is because at higher area-fractions the particle phase produces stronger hydrodynamic resistance for the flow of the continuous-phase fluid in the space between the particles. (b) (a) 𝜒 𝑐 = 1 + (1 − 𝜙̅ 𝑎 ) 𝑓̅ ( -9) where 𝑓 ̅ = 𝑓 𝜇 𝑠 /ℎ is a dimensionless friction factor. An increase in the area fraction 𝜙 ̅ 𝑎 leads to a decrease in the dimensionless friction factor 𝑓 ̅ as well as its prefactor 𝑎 , causing 𝜒 𝑐 to grow towards unity. Fig. 5a shows the values of 𝜈̅ 𝑝𝑝 versus the area fraction, evaluated at locations 1 and 2. As expected, the suspension permeability coefficient 𝜈̅ 𝑝𝑝 is a decaying function of the area fraction 𝜙̅ 𝑎 , because of the increased flow resistance at higher droplet concentrations. The dashed line in Fig. 5a is the Clausius-Mossotti approximation for permeability coefficients 𝜈̅ 𝑝𝑝 𝜈̅ CM𝑝𝑝 = 𝜈̅ (1 − 0.5 𝜇̃ 𝜙 𝑎 )(1 + 0.5 𝜇̃ 𝜙 𝑎 ) ( -10) where the low-density expansion coefficient 𝜇̃ = 3.7 is estimated from the experiment with the lowest area fraction (details in the ESI Section 4). The droplet mobility coefficient 𝜈̅ 𝑡𝑝 evaluated from the measured transport coefficients 𝜒 𝑐 and 𝜈̅ 𝑝𝑝 using eqn (3-4) is depicted in Fig. 5b. The coefficient 𝜈̅ 𝑡𝑝 has a weaker dependence on the area fraction compared to 𝜈̅ 𝑝𝑝 , because the convective transport coefficient 𝜒 𝑐 is an increasing function of area fraction, partially compensating for the decay of 𝜈̅ 𝑝𝑝 in eqn (3-4). During the capillary imbibition of disk-like droplets, a sharp front forms between the suspension domain and droplet-free region (Figs. 1b (i) and 2). We also observe that the density fluctuations in the suspension exhibit wave-like dynamics (Fig. 6 and ESI video). Both the density-front formation and the propagating density waves can be explained using the convective transport eqn (3-3) and our measurements showing that the transport coefficient 𝝌 𝒄 is an increasing function of 𝝓̅ 𝒂 , Fig. 4. Domains with higher area fraction move faster than lower-density regions. In a zone where the area fraction decreases, such as the leading end of the suspension region, the higher density portion of the suspension catches up with the lower density portion. This results in sharpening the area-fraction gradient and eventually forming a density front. For moderate amplitude density fluctuations, a similar mechanism produces density waves propagating in the forward direction, because the leading end of the dense region collects particles.
Quantitatively, front formation and wave propagation are described by the convective-transport equation for the suspension density. This equation, obtained by combining the constitutive relation eqn (3-3) with the continuity condition for the particle phase, 𝜕𝜙̅ 𝑎 𝜕𝑡 = −∇ || (𝜙̅ 𝑎 𝑈̅) ( -11) can be written in the form of the generalized inviscid Burgers equation, 𝜕𝜙̅ 𝑎 𝜕𝑡 = −𝑐(𝜙̅ 𝑎 )∇ || 𝜙̅ 𝑎 ( -12) where 𝑐(𝜙̅ 𝑎 ) = [𝝌 𝒄 (𝜙̅ 𝑎 ) + 𝜙̅ 𝑎 d𝝌 𝒄 (𝜙̅ 𝑎 )d𝜙̅ 𝑎 ] 𝑢̅ ( -13) is the area-fraction-dependent velocity of the density waves. The Burgers equation eqn (3-12) predicts the formation of density shocks between the high- and low-density regions for systems where the wave velocity 𝒄(𝝓̅ 𝒂 ) is an increasing function of the area fraction 𝝓̅ 𝒂 . This is the case for our experimental system. Shock formation can be demonstrated using the method of characteristics, and it results from the characteristics crossing. The observed sharp front between the suspension and clear-fluid regions is an example of a density shock predicted by the Burgers equation eqn (3-12). Fig. 5.
Transport coefficients of disk-like droplets. (a) The effective permeability coeficinets 𝜈̅ 𝑝𝑝 of suspension of disk-like droplets vs 𝜙̅ 𝑎 . Square markers are for location 1 and diamond markers are for location 2. The dashed line is Clausius-Mossoti approximation where expansion coefficient 𝜇̃ = 3.7 is estimated from the experiment with the lowest area fraction. 𝜈̅ 𝑝𝑝 is a decaying function of 𝜙̅ 𝑎 . (b) The droplet mobility coefficient 𝜈̅ 𝑡𝑝 vs 𝜙̅ 𝑎 . The coefficient 𝜈̅ 𝑡𝑝 has a weaker dependence on the area fraction compared to 𝜈̅ 𝑝𝑝 , because the convective transport coefficient 𝜒 𝑐 is an increasing function of area fraction, partially compensating for the decay of ν ̅ ^pp in eqn (3-4). Is (b)(a) Fig. 6.
Density wave propagation during the imbibition of disk-like droplets.
The wave velocity predicted by eqn (3-13) is greater than average droplet velocity. Therefore, the density waves associated with spontaneously occurring density fluctuations propagate faster than the average particle velocity. As a result, the distance between the front of the density wave (depicted by dashed line) and the front of the suspension region (solid vertical blue line) is decreasing. e note that the wave velocity eqn (3-13) is greater than the convective velocity eqn (3-3). Therefore, the density waves associated with spontaneously occurring density fluctuations propagate faster than the average particle velocity and t the suspension front. This behavior is shown in Fig. 6, where the distance between the propagating low-density region and the front gradually decreases. The observed wave-propagation direction relative to the average particle velocity in the quasi-2D imbibition system is opposite the one reported for 1-D drop trains . For linear chains oriented parallel to the flow direction, the particle mobility decreases with the increasing particle density per unit length. The decrease of the mobility coefficient explains the backward direction of the propagating density wave. For pressure driven flow of disk-like droplets in a microfluidic Hele – Shaw cell, density-wave propagation and shock formation was reported by Beatus, Tlusty, and Bar – Ziv. Our analysis provides a generalization of their results in the context of the imbibition process.
As shown in Figs. 7 and 8, droplets self-organize into a variety of transient local arrangements, which include short chains of different orientations, some elongated and some more compact, as well as high density regions separated by gaps of much lower particle density. Our numerical simulations using the dipolar-interactions model described in section 2.7 capture the emergence of these transient structures, as depicted in Fig. 7. The simulations also replicate the sharp border between the suspension and clear fluid regions. Apart from the axisymmetric near-field repulsion, particles in our model interact only via the Hele – Shaw dipoles. The agreement between the experimental and simulated microstructural features indicates that the dipolar hydrodynamic interactions provide the main organizing factor that leads to the emergence of the observed microstructural dynamics. The chains generated during the imbibition process are usually short, with a typical length of several particles. They have relatively short lifetimes because of inherent instabilities or collisions with other microstructural features. We note that isolated transverse chains are more stable than longitudinal ones . In general, the microstructure generated during the imbibition is much less ordered than the one seen for a system of share-driven deformable dropletss for which dynamics are governed by the Hele-Shaw quadrupolar interactions instead of dipoles . As illustrated in the time-lapse sequence of video images presented in Fig. 8a, chains and other elongated clusters with transverse orientation (across the channel) move faster than longitudinally oriented ones. Fig. 8b shows that diagonally oriented clusters have a nonzero transverse velocity component. Both the variation of the chain/cluster velocity with its orientation and the higher velocity of dense regions discussed in Sec. 3.2.1 contribute to the large particle velocity fluctuations that occur during the imbibition process, Fig. 9. The cluster velocity variation and its associated cluster collisions, Fig. 7a, also give rise to high density regions, enhancing the area-fraction fluctuations. The Hele – Shaw dipolar model captures these features of droplet dynamics. The variation of the velocity and direction of cluster motion stems from the geometry of dipolar backflow patterns, as depicted in the schematic shown in Fig. 10, panels (ii). To further elucidate these phenomena, we approximate a particle chain as an elongated rigid body. In this simple approximation, chain motion is described by the mobility relation
𝑽 = 𝝂 𝑐 · 𝒖 ( -14) where 𝒖 is the velocity of the incident flow, Fig 7. Suspension model with dipolar Hele-Shaw hydrodynamic interactions. (a) Similar microstructural features between simulation and experiment. The meniscus is drawn as a guide for the eye (b) Comparison between predicted microstructural features (bottom images) and experiments (top images). The same geometrical shape is used to highlight similar structural elements in experimental snapshots and simulations.
Fig 8. Dynamics of longitudinal, transverse, and diagonally oriented droplet chains. (a) The collision of a fast moving transverse droplet chain with a slow moving longitudinal chain resulting in a scattering of the droplets in (i) imbibition experiment, (ii) simulations based on dipolar interactions (b) The transverse motion of diagonal droplet chains toward the side walls in (i) capillary imbibition experiments, (ii) simulations based on dipolar interactions. 𝑐 = 𝜈 ||𝑐 𝒕̂𝒕̂ + 𝜈 ⊥𝑐 𝒏̂𝒏̂ ( -15) is the chain mobility tensor, 𝒕̂ and 𝒏̂ are unit the vectors parallel and perpendicular to the chain, and 𝜈 ||𝑐 and 𝜈 ⊥𝑐 are the longitudinal and transverse mobility components. In strongly confined systems, the mobility tensor of a particle chain or an elongated particle cluster is strongly anisotropic
34, 50 , with the transverse mobility significantly larger than the longitudinal mobility, 𝜈 ⊥𝑐 > 𝜈 ||𝑐 ( -16) In a quasi-2D system, the fluid cannot easily bypass a chain oriented normal to the flow direction, thus pushing the chain forward. In contrast, a longitudinal chain can be more easily bypassed by the flow, Fig. 10 panels (i). Combining relations (3-14) and (3-15) we get the expressions 𝑉 𝑥 = [𝜈 ||𝑐 sin (𝜃) + 𝜈 ⊥𝑐 cos (𝜃)]𝑢̅ ( -17a) and 𝑉 𝑦 = (𝜈 ⊥𝑐 − 𝜈 ||𝑐 )sin(2𝜃)𝑢̅ ( -17b) for the longitudinal and transverse chain velocity components 𝑉 𝑥 and 𝑉 𝑦 , where 𝜃 is the angle between the flow direction 𝒆̂ 𝑥 and the normal to the chain 𝒏̂ . Taking into account the anisotropy of the chain mobility tensor, eqn (3-16), relation (3-17a) shows that chains with approximately transverse orientation 𝜃 ≈ 0 ° move faster along the channel than the longitudinal ones ( 𝜃 ≈ 90 °). Moreover, diagonally oriented chains (0°< 𝜃 <90°) have a nonzero transverse velocity component, predicted by eqn (3-17b). The experimentally observed chain behaviour is consistent with these theoretical results. Imbibition of spherical droplets 3.3.1 Droplet distribution and overall imbibition dynamics
Imbibition of a monolayer of spherical droplets differs considerably from imbibition of strongly confined disk like droplets. The difference stems primarily from the opposite relation between the droplet-phase and suspension velocities 𝑈̅ and 𝑢̅ . Spherical droplets move faster than the suspension (i.e., we have 𝜒 𝑐 > 1 ), whereas disk like droplets move slower ( 𝜒 𝑐 <1 ), as shown in Sec. 3.2. Moderately confined spherical droplets are faster than the average flow, because they reside in the midplane of the channel and therefore sample the streamlines with higher velocities. The focusing of droplets has been attributed to inertial migration or deformation-based migration
51, 52 . As a result of the faster motion, spherical droplets reach the meniscus and form a dense droplet-band, as seen in the bottom panel of Fig. 1b and in Fig. 11. In contrast, there is a clear-fluid region behind the meniscus for the disk-like droplets. The dense
Fig 10. Schematics showing the flow of fluid streamline and dipolar backflows in longitudinal, transverse and diagonal droplet chains. (i) Fluid streamlines around (a) longitudinal, (b) transverse, and (c) diagonally oriented chains. (ii)
Geometry of a dipolar backflow pattern for a) longitudinal (b) transverse (c) diagonally oriented chain.
Fig. 9 – Distribution of x and y component of the droplet velocity
The distribution of the velocity of the droplets in the x-direction 𝑣 𝑥 (the flow direction) is wider than the distribution in the y-direction ( 𝑣 𝑦 ) (parallel to flow direction. This is due to the variation in the speed of droplet chain with different orientations as well as wave propagations discussed in Section 2.3.2.1. Fig. 11 – Distribution of spherical droplets in the channel . The area fraction of emulsion is measured in each of 13 equal divisions of a portion of the channel. The bottom panel shows the measured values of area fraction. We observe an axial gradient in the area fraction of the emulsion as well as large fluctuations in the area fraction. The snapshot is from an experiment with average area fraction for an emulsion with average area fraction of 𝜙̅ 𝑎 =0.53. egion in the spherical drop system often has close-packed hexagonal domains and disordered regions. The border between these regions is unstable and diffuse. There is a systematic concentration droplet away from the meniscus, as quantitatively shown in the bottom panel of Fig. 11. This gradually varying particle density distribution qualitatively differs from the one observed for disk-like drops, for which there is a sharp density front and no systematic density gradient. Moreover, according to the bottom panels of Figs. 2 and 11 the area fraction fluctuations are considerably larger for spherical droplets than those that occur for disk-like drops. Fig. 12 shows the imbibition time for spherical droplets at locations 1 and 2 of the channel. The imbibition time is an increasing function of the area fraction, although this dependence is weaker than the one observed for disk-like droplets. There are also larger fluctuations from experiment to experiment. The fluctuations result from large emulsion density fluctuations and from the presence of the dense droplet region next to the meniscus, which affects the contact angle and the capillary pressure driving the imbibition. The density fluctuations and variation of the contact angle also cause inaccuracies in estimating the time for the meniscus to reach location 2; hence the difference of the imbibition time evaluated using observations at different locations. The dynamics of the suspension of spherical droplets is governed by the same set of constitutive equations, (3-1) to (3-3), as those that describe disk-like droplets. In the disk-like droplet imbibition, the presence of a uniform suspension region and a clear-fluid region separated by a sharp front allow determining the transport coefficients 𝝌 𝒄 , 𝝂̅ 𝒑𝒑 , and 𝝂̅ 𝒕𝒑 with good accuracy. For spherical droplets, however, significant area-fraction nonuniformities and the lack of a front between the suspension and the clear fluid limits the ability to evaluate the transport coefficients. Therefore, the results are less complete than those given in Sec. 3.2. We first consider the convective transport coefficient 𝜒 𝑐 . Its evaluation requires measurement of the volumetric suspension velocity 𝑢̅ and the particle-phase velocity 𝑈̅ , eqn (3-3). As in the case of disk-like droplets, the suspension velocity 𝑢̅ can be determined from the motion of the meniscus. However, the particle-phase velocity 𝑈̅ is more difficult to measure because there is no suspension front moving with the velocity of the particle phase. In the low-area fraction experiments, 𝑈̅ can be determined by particle tracking, but at higher area fractions, the droplets in high-density cluster are not easily tracked. Therefore, we provide results only for the lowest area fraction system, using data from experiments imaged at location 2. Measurements yield 𝜒 𝑐 = 1.26 for 𝜙̅ 𝑎 = 0.2 and 𝜒 𝑐 = 1.29 for, 𝜙̅ 𝑎 = 0.23 . As expected, 𝜒 𝑐 > 1 , consistent with the qualitative observations. The experimentally determined value is lower by approximately 15% than the result 𝜒 =1.5 for a point particle the midplane of the channel moving in the parabolic flow 𝑢(𝑧) = 16 𝑧ℎ (𝑧ℎ − 1) 𝑢̅ ( -18) where the channel walls are at 𝑧 = 0 and 𝑧 = ℎ . The measured values are also lower, by less than 10%, than the value 𝜒 𝑐Faxen =1.39 that is obtained by applying Faxen’s law 𝑈̅ = 𝑢 + 𝑑
24 ∇ 𝑢 ( -19) to an isolated solid sphere of diameter 𝑑 = 0.65ℎ moving in the parabolic flow, eqn (3-18). The drop velocity cannot exceed the point-particle velocity, and it is also unlikely to be larger than the Faxen’s prediction , in which frictional forces due to the wall presence are neglected. We conclude that the coefficient 𝜒 𝑐 ( 𝜙̅ 𝑎 ) remains nearly constant or decreases when the area fraction 𝜙̅ 𝑎 is increased. This conjecture is consistent with the theoretical calculations for a hard-sphere monolayer, which has very weak area fraction dependence of 𝜒 𝑐 (Sec. 4.2). It follows that the suspension of spherical drops does not show shocks and density-wave propagation that are observed for disk-like droplets. Evaluation of the transport coefficient 𝜈̅ 𝑝𝑝 requires measuring the volumetric suspension velocity 𝑢̅ and determining the local pressure gradient in a uniform suspension region. While 𝑢̅ can be obtained from the meniscus velocity, the local pressure gradient cannot be evaluated from our experimental data because of suspension nonuniformities. Therefore, we calculate the average value of 𝜈̅ 𝑝𝑝 for the entire channel: 𝑢̅ = 𝜈̅ avg𝑝𝑝 ∇ || P avg ( -20) where ∇ || 𝑃 avg = Δ𝑃/𝐿 is the average pressure gradient, with the pressure drop Δ𝑃 evaluated from the Young-Laplace equation. This is similar to treatment of disk-like droplets. The average permeability coefficient 𝜈̅ avg𝑝𝑝 , based on measurements performed at locations 1 and 2, is plotted in Fig. 13 versus the average area fraction. The data show that 𝜈̅ 𝑎𝑣𝑔𝑝𝑝 is a decaying function of the area fraction. As discussed in Sec. 4.2, Fig. 12.
Filling time of capillary tube for spherical droplets.
The time required to fill the capillary tube with spherical emulsion, measured at (a) Location 1 and (b) Location 2 of the channel. The filling time increases with area fraction because the total friction force opposing the imbibition is larger at higher area fractions. The solid lines are least-square fits drawn as a guide for the eyes. (b)(a) he decay is much stronger than the one expected according to a comparison with results for a hard-sphere monolayer. Based on an analysis of high magnification images (Fig. 14), we attribute this discrepancy to displacements of highly packed drops out of the midplane of the channel. The displaced drops can create significant friction with the walls, hindering suspension flow. Our average results thus imply a high permeability contrast between low-density and nearly close-packed domains of the channel. By similar arguments, we also expect that the coefficient 𝜒 𝑐 may significantly decrease at high area fractions. Strong fluctuations of the area fraction of spherical droplets and the existence of a systematic area-fraction gradient in the axial direction, stem from the instability of the dense droplet-band, which is locally displaced by less-dense suspension. An example of this instability is presented in the time-lapse image sequence depicted in Fig. 15. The observed instability is analogous to the Saffman – Taylor viscous fingering instability in a Hele – Shaw cell. The fingering occurs when a high viscosity fluid is displaced by a fluid with lower viscosity. Since the parts of the channel filled with the low-viscosity fluid are more permeable, fluid flows more easily in these domains, creating low-viscosity fingers. A similar behavior occurs here, where low droplet area-fraction regions have a larger channel permeability coefficient 𝜈̅ 𝑝𝑝 than the high area-fraction regions. See the average results shown in Fig. 13. Hydrodynamic stresses associated with the motion of low-density fingers can result in meniscus deformation, Fig. 15. The deformation in turn facilitates finger formation because it breaks the hexagonally ordered suspension structure in the closely packed region near the suspension front. Closely packed ordered emulsion has a nonzero yield stress, so development of an ordered structure may prevent finger formation, allowing a buildup of a large dense domain. This behavior is seen in Fig. S7 Fig. 13. Average suspension permeability coeficient 𝝂̅ 𝒂𝒗𝒈𝒑𝒑 for spherical droplets. The measured values of 𝜈̅ 𝑎𝑣𝑔𝑝𝑝 for disk-like droplets at (a) location 1, and (b) location 2. 𝜈̅ 𝑎𝑣𝑔𝑝𝑝 is a decaying funcion of area fraction. However, the decay is much stronger compared to numerical sumulations for hard-sphere monolayers (Section 4.2). (b)(a) Fig. 14. High magnification images of dense regions formed during the imbibition of spherical droplets.
High magnification images (4X) of dense regions formed during the imbibition of spherical droplets at (a) upstream of the meniscus (b) at the meniscus. The regions marked by blue circles show droplets with overlapping boundaries, indicating the displacement of the droplets out of midplane.
Fig. 15. Fingering instability and interface deformations during the imbibition of spherical droplets.
Time-lapsed images from an experiment with 𝜙̅ 𝒂 = (i)(ii) (iii) (iv) t = 124 ms t = 0 mst = 281 ms t = 680 ms Fig. 16. Hele-Shaw dipolar model predicts the formation of hexagonally packed droplet-band.
Qualitative comparison between the droplet band formed in an experimnt with 𝜙̅ 𝒂 = n the ESI, which shows results at a higher confinement ratio. Formation of hexagonal structures is also predicted by our numerical Hele – Shaw dipolar model with meniscus approximated using a repulsive force acting in the direction opposite to the flow, Fig. 16. When a low-density finger approaches the meniscus it usually turns towards one of the side walls of the channel. After touching the wall, the finger causes separation of a dense droplet cluster from the close-packed droplet band near the meniscus. As seen from the time-lapse images presented in Fig. 11, the separated clusters migrate away from the meniscus, gradually decaying from random droplet motions. This process periodically repeats, causing the observed large area- fraction fluctuations and formation of a systematic area-fraction gradient seen in Fig. 11. Fingering instabilities were reported in pressure driven flow of a suspension in a radial Hele-Shaw cell
36, 37 . Our results demonstrate that the instability also occurs during imbibition and that it can significantly affect the imbibition dynamics. Discussion and conclusions
The results presented in Sec. 3.2 and 3.3 reveal significant differences between the imbibition hydrodynamics in the disk-like- and spherical-droplet systems. The observed differences are most evident near the meniscus, where either clear-fluid or densely-packed region forms. However, we also found fundamental differences in the mesoscale droplet-distribution dynamics away from the meniscus.
For disk-like droplets, there is a large variation between the velocities of longitudinally and transversely oriented chains. Such variation is less pronounced for spherical droplets, consistent with the fact that the mobility tensor (eqn (3-15)) is less anisotropic for weaker confinement
33, 34 . We also note that transverse chains of spherical droplets are less stable than the transverse chains/clusters of strongly confined disk-like droplets, as predicted by the analysis of chain stability at different confinement . Dense clusters of disk-like droplets often span the width of the channel and exhibit density wave dynamics, as discussed in Sec. 3.2.3. Both the droplets in these clusters and the waves move faster than the drops in the low-density regions. In contrast, for spherical droplets, dense clusters usually span only part of the channel width and move slower than the low-density regions. Moreover, no density waves are observed in the spherical-droplet system. The origin of distinct dynamics of density fluctuations is explained in the schematic shown in Fig. 17. Due to the flow continuity condition, the dynamics of transversely oriented clusters of disk-like-droplets (Fig. 17a) are primarily controlled by the convective trasport equation (3-3). The drop density varies mainly in the longitudinal direction. Because of the continuity condition, the average suspension velocity is the same in the region behind a cluster, in the cluster, and in front of the cluster. The convective transport coefficient has the largest value in the dense region, 𝜒 𝑐1 , 𝜒 𝑐3 < 𝜒 𝑐2 , so one concludes that the dense cluster moves faster than the lower area-fraction regions ( 𝑈̅ , 𝑈̅ < 𝑈̅ ). In a system of spherical droplets, the clusters that separate from the dense droplet-band do not span the entire width of the channel. As explained in Fig. 1 b, their motion relative to the neighboring low-density regions is thus governed by the mobility relations (3-1) and (3.2). Since the pressure gradient in the cluster and a latelarly adjacent low-density domain are approxiamtely the same, and the dense region has a lower permeability coefficient 𝜈̅ 𝑎𝑣𝑔,1𝑝𝑝 < 𝜈̅ 𝑎𝑣𝑔,2𝑝𝑝 , it moves slower than the low-density region ( 𝑢̅ < 𝑢̅ ). Noting that the convective transport coefficient 𝜒 𝑐 is reduced in dense clusters due to out-of-plane particle displacements, the particle velocity is further diminished. As a consequence, a significant motion of clusters away from meniscus is generated. This backward motion results in particle transport away from the meniscus, balancing the forward motion of droplets in low density regions, and preventing a continuous buildup of the densely packed suspension domain behind the meniscus. In this section we compare the experimentally determined transport coefficients 𝜒 𝑐 , 𝜈̅ 𝑝𝑝 , and 𝜈̅ 𝑡𝑝 to our benchmark calculations for a confined hard-sphere monolayer driven by a pressure gradient in a parallel-wall channel. The benchmark calculations were performed using the Cartesian representation method for a hexagonal lattice of spherical particles moving in the midplane of the channel. The details of our technique along with a theoretical analysis of particle transport in parallel-wall channels are provided in . The calculations are Fig. 17. Comparing the dynamics of dense droplet clusters for disk-like and spherical droplet systems. (a) The dense clusters of disk-like droplets usually span the whole width of the channel and its dynamics is governed by 𝜒 𝑐 resulting in faster motion of the dense region (b) The dense cluster of spherical droplets do not span the whole width of channel and their dynamics is governed by 𝜈̅ 𝒂𝒗𝒈𝒑𝒑 resulting in a slower motion of the dense region. Disk-like dropletsSpherical droplets a (b) erformed for a hexagonal lattice, while the experimental systems are mostly disordered. However, based on the comparison between hexagonal-lattice and square-lattice results, we expect that the microscale particle distribution does not strongly affect transport coefficients, as long as there are no large area fraction fluctuations (which in our calculations are not present). The theoretical hard-sphere results for 𝜒 𝑐 , 𝜈̅ 𝑝𝑝 , and 𝜈̅ 𝑡𝑝 at confinement rations ranging from 𝑑/ℎ = 0.33 to 𝑑/ℎ = 0.9999 are depicted in Fig. 18 (lines) along with the corresponding experimental data for disk-like and spherical droplet systems (symbols). The results are plotted vs the area fraction 𝜙̅ 𝑎 , and each line represents calculation for a different confinement ratio. We find that the measured transport coefficients for disk-like droplets quantitatively agree with the theoretical hard-sphere calculations for a strongly confined particle monolayer. Specifically, the measurements performed at locations 1 and 2 can be mapped onto the hard-sphere results for ℎ/𝑑 = 0.999 and ℎ/𝑑 = 0.9999 , respectively. The spherical-droplet system overall agreement is less satisfactory, except at the lowest area fraction 𝜙̅ 𝑎 ≈ 𝜒 𝑐 agrees well with the hard-sphere results for the same confined ratio as the one used in our experiments. The overall agreement between experiment and benchmark calculations indicates particle motion in both the hard-sphere and droplet systems is governed by similar hydrodynamic mechanisms. As argued in Sec. 3.2.2 and the ESI, section 5, particle dynamics at strong confinements is controlled by the competition between the pressure force pushing the particles along the channel and the particle – wall hydrodynamic friction. Moreover, the fluid motion relative to the particle phase is governed by the inter-phase permeability coefficient 𝜇 𝑠 , which depends on particle spacing. For disk-like droplets, the particle-wall friction is generated in the thin fluid film between the wall and the drop interface, and for highly confined rigid spheres, it is produced in the lubrication layer in the near-contact region. Since the hard-sphere lubrication region is much smaller than the flattened-droplet film, we need to use very strong confinements to match the theoretical results to the experimental ones. The confinement ratio ℎ/𝑑 required to match the results obtained from measurements at location 1 is smaller than the confinement ratio needed for location 2. This is consistent with our conjecture that the location-dependence of the measured values of the transport coefficients stems from the increased film thickness at the beginning of the imbibition process due to the faster drop motion. Similar area-fraction dependence of the droplet and hard-sphere transport coefficients seen in Fig. 18 suggests that the backflow associated with dipolar Hele-Shaw interactions is a key factor in suspension dynamics for both systems. For spherical droplets in the moderately confined system, we find a good agreement between the experimental and theoretical results only at low area fractions. The agreement is especially good for the convective transport coefficient 𝜒 𝑐 , for which accurate particle-tracking results are available for the area fraction 𝜙̅ 𝑎 ≈ 𝜈̅ avg𝑝𝑝 , estimated by averaging over the entire non-uniform suspension region. We find that the theoretical and experimental results significantly differ at higher values of fraction 𝜙̅ 𝑎 . According to our theoretical results, the coefficient 𝜈̅ 𝑝𝑝 decreases by less than 10% in area-fraction domain for which experimental results are available, whereas the experimentally determined values of 𝜈̅ avg𝑝𝑝 , decrease by almost 40%. This significant difference between the experimental and theoretical results stems, most likely, from the displacements of spherical drops out of the midplane of the channel in high area-fraction regions, as discussed in Sec. 3.3.2. The presented results demonstrate rich phenomenology of the emulsion imbibition process under moderate and strong confinements. In particular we have shown that the drop distribution can be controlled by the degree of confinement. Under strong confinement where the droplets deform to a disk-like shape, we observe formation of droplet-free region behind the meniscus. Moreover, the suspension region has a sharp front boundary due to the presence of shockwaves. Droplets in the suspension region show strong velocity and arae fraction fluctuations with density wave propagations. During the imbibition of moderately confined spherical droplets, a dense droplet-band forms behind the meniscus that undergoes fingering instability. As opposed to disk-like droplets there is no sharp boundary between droplet-band and suspension
Fig. 18. Experimental transport coefficients 𝝌 𝒄 , 𝝂̅ 𝒑𝒑 , and 𝝂̅ 𝒕𝒑 compared with numerical simulations for confined hard spheres. Lines represent the values from numerical simulations based-on confined hard-spheres and symbols represent the experimental measurements for (a) 𝜒 𝑐 , (b) 𝜈̅ 𝑝𝑝 , and (c) 𝜈̅ 𝑡𝑝 . For disk-like droplets, the experimental measurements match well with the hard-sphere results for ℎ/𝑑 = 0.999 and ℎ/𝑑 = 0.9999 . For spherical-droplets measurements agree only at lower area fractions . (a) (b) (c) egion. However, our experiments with spherical droplets at stronger confinements shows formation of closely-packed dense regions that can prevent fingering instability. Acknowledgements
Acknowledgement is made to the donors of the American Chemical Society Petroleum Research Fund for support of M.W.V and M.N in this research. S.S. and J.B. were supported by NSF Grant CBET 1603627. We thank Purushottam Soni for the surface tension measurement.
References
1. Delamarche, E.; Bernard, A.; Schmid, H.; Michel, B.; Biebuyck, H., Patterned delivery of immunoglobulins to surfaces using microfluidic networks.
Science (5313), 779-781. 2. Gökçe, O.; Castonguay, S.; Temiz, Y.; Gervais, T.; Delamarche, E., Self-coalescing flows in microfluidics for pulse-shaped delivery of reagents.
Nature (7777), 228-232. 3. Morrow, N. R.; Mason, G., Recovery of oil by spontaneous imbibition.
Current Opinion in Colloid & Interface Science (4), 321-337. 4. Clarke, A.; Blake, T.; Carruthers, K.; Woodward, A., Spreading and imbibition of liquid droplets on porous surfaces. Langmuir (8), 2980-2984. 5. Libri, G., Histoire des sciences mathématiques en Italie: depuis la renaissance des lettres jusqu'à la fin du dix-septieme siècle . HW Schmidt: 1865; Vol. 2. 6. Washburn, E. W., The dynamics of capillary flow.
Physical review (3), 273. 7. Quéré, D., Inertial capillarity. EPL (Europhysics Letters) (5), 533. 8. Fries, N.; Dreyer, M., The transition from inertial to viscous flow in capillary rise. Journal of colloid and interface science (1), 125-128. 9. Bracke, M.; De Voeght, F.; Joos, P., The kinetics of wetting: the dynamic contact angle. In
Trends in Colloid and Interface Science III , Springer: 1989; pp 142-149. 10. Siebold, A.; Nardin, M.; Schultz, J.; Walliser, A.; Oppliger, M., Effect of dynamic contact angle on capillary rise phenomena.
Colloids and surfaces A: Physicochemical and engineering aspects (1), 81-87. 11. Popescu, M. N.; Ralston, J.; Sedev, R., Capillary rise with velocity-dependent dynamic contact angle.
Langmuir (21), 12710-12716. 12. Zhou, R.; Chang, H.-C., Capillary penetration failure of blood suspensions. Journal of colloid and interface science (2), 647-656. 13. Chakraborty, S., Dynamics of capillary flow of blood into a microfluidic channel.
Lab on a Chip (4), 421-430. 14. Berli, C. L.; Urteaga, R., Asymmetric capillary filling of non-Newtonian power law fluids. Microfluidics and nanofluidics (6), 1079-1084. 15. Yao, Y.; Alexandris, S.; Henrich, F.; Auernhammer, G.; Steinhart, M.; Butt, H.-J.; Floudas, G., Complex dynamics of capillary imbibition of poly (ethylene oxide) melts in nanoporous alumina. The Journal of Chemical Physics (20), 203320. 16. Beard, J. D.; Rouholamin, D.; Farmer, B. L.; Evans, K. E.; Ghita, O. R., Control and modelling of capillary flow of epoxy resin in aligned carbon nanotube forests.
RSC Advances (49), 39433-39441. 17. Psaltis, D.; Quake, S. R.; Yang, C., Developing optofluidic technology through the fusion of microfluidics and optics. Nature (7101), 381-386. 18. Dhar, J.; Jaggi, P.; Chakraborty, S., Oscillatory regimes of capillary imbibition of viscoelastic fluids through concentric annulus.
RSC advances (65), 60117-60125. 19. Sumanasekara, U. R.; Azese, M. N.; Bhattacharya, S., Transient penetration of a viscoelastic fluid in a narrow capillary channel. Journal of Fluid Mechanics , 528-552. 20. Shou, D.; Fan, J., The fastest capillary penetration of power-law fluids.
Chemical Engineering Science , 583-589. 21. Siddique, J.; Anderson, D., Capillary rise of a non-Newtonian liquid into a deformable porous material.
Journal of Porous Media (12). 22. Srivastava, N.; Burns, M. A., Analysis of non-Newtonian liquids using a microfluidic capillary viscometer. Analytical chemistry (5), 1690-1696. 23. Rashid, M. A. I.; Ansari, S.; Waghmare, P. R.; Nobes, D., Capillary filling of Non-Newtonian fluids. 24. Holloway, W.; Aristoff, J. M.; Stone, H. A., Imbibition of concentrated suspensions in capillaries. Physics of Fluids (8), 081701. 5. Beatus, T.; Bar-Ziv, R. H.; Tlusty, T., The physics of 2D microfluidic droplet ensembles. Physics reports (3), 103-145. 26. Beatus, T.; Tlusty, T.; Bar-Ziv, R., Burgers shock waves and sound in a 2d microfluidic droplets ensemble.
Physical review letters (11), 114502. 27. Singha, S.; Malipeddi, A. R.; Zurita-Gotor, M.; Sarkar, K.; Shen, K.; Loewenberg, M.; Migler, K. B.; Blawzdziewicz, J., Mechanisms of spontaneous chain formation and subsequent microstructural evolution in shear-driven strongly confined drop monolayers.
Soft matter (24), 4873-4889. 28. Saffman, P., The lift on a small sphere in a slow shear flow. Journal of fluid mechanics (2), 385-400. 29. Baron, M.; Bławzdziewicz, J.; Wajnryb, E.,
Hydrodynamic crystals: Collective dynamics of regular arrays of spherical particles in a parallel-wall channel.
Physical review letters (17), 174502. 30. Beatus, T.; Tlusty, T.; Bar-Ziv, R., Phonons in a one-dimensional microfluidic crystal.
Nature Physics (11), 743. 31. Hashimoto, M.; Garstecki, P.; Stone, H. A.; Whitesides, G. M., Interfacial instabilities in a microfluidic Hele-Shaw cell. Soft Matter (7), 1403-1413. 32. Cui, B.; Diamant, H.; Lin, B.; Rice, S. A., Anomalous hydrodynamic interaction in a quasi-two-dimensional suspension. Physical review letters (25), 258301. 33. Bhattacharya, S.; Bławzdziewicz, J.; Wajnryb, E.,
Many-particle hydrodynamic interactions in parallel-wall geometry: Cartesian-representation method.
Physica A: Statistical Mechanics and its Applications (2-4), 294-340. 34. Bhattacharya, S.; Blawzdziewicz, J.; Wajnryb, E., Hydrodynamic interactions of spherical particles in suspensions confined between two planar walls. arXiv preprint physics/0505056 . 35. Janssen, P.; Baron, M.; Anderson, P.; Blawzdziewicz, J.; Loewenberg, M.; Wajnryb, E., Collective dynamics of confined rigid spheres and deformable drops.
Soft Matter (28), 7495-7506. 36. Tang, H.; Grivas, W.; Homentcovschi, D.; Geer, J.; Singler, T., Stability considerations associated with the meniscoid particle band at advancing interfaces in Hele-Shaw suspension flows. Physical Review Letters (10), 2112. 37. Kim, J.; Xu, F.; Lee, S., Formation and destabilization of the particle band on the fluid-fluid interface. Physical review letters (7), 074501. 38. Anna, S. L.; Bontoux, N.; Stone, H. A., Formation of dispersions using “flow focusing” in microchannels.
Applied physics letters (3), 364-366. 39. Dreyfus, R.; Tabeling, P.; Willaime, H., Ordered and disordered patterns in two-phase flows in microchannels. Physical review letters (14), 144505. 40. Xia, Y.; Whitesides, G. M., Soft lithography. Annual review of materials science (1), 153-184. 41. Fuguet, E.; Ràfols, C.; Rosés, M.; Bosch, E., Critical micelle concentration of surfactants in aqueous buffered and unbuffered systems. Analytica Chimica Acta (1-2), 95-100. 42. Rasband, W. S., ImageJ. Bethesda, MD: 1997. 43. Kalman, R. E., A new approach to linear filtering and prediction problems.
Journal of basic Engineering (1), 35-45. 44. Kuhn, H. W., The Hungarian method for the assignment problem. Naval research logistics quarterly -97. 45. Bławzdziewicz, J.; Wajnryb, E., Equilibrium and nonequilibrium thermodynamics of particle-stabilized thin liquid films.
The Journal of chemical physics (19), 194509. 46. Beatus, T.; Bar-Ziv, R.; Tlusty, T., Anomalous microfluidic phonons induced by the interplay of hydrodynamic screening and incompressibility.
Physical review letters (12), 124502. 47. Bławzdziewicz, J.; Wajnryb, E., An analysis of t he far-field response to external forcing of a suspension in the Stokes flow in a parallel-wall channel.
Physics of Fluids (9), 093303. 48. Blawzdziewicz, J.; Goodman, R.; Khurana, N.; Wajnryb, E.; Young, Y.-N., Nonlinear hydrodynamic phenomena in stokes flow regime. Physica D: Nonlinear Phenomena (14), 1214-1224. 49. Migler, K. B., String formation in sheared polymer blends: Coalescence, breakup, and finite size effects.
Physical review letters (6), 1023. 50. Bhattacharya, S. ; Bławzdziewicz, J.; Wajnryb, E., Hydrodynamic interactions of spherical particles in Poiseuille flow between two parallel walls.
Physics of Fluids (5), 053301. 1. Olla, P., Simplified model for red cell dynamics in small blood vessels. Physical review letters (2), 453. 52. Goldsmith, H.; Mason, S., The flow of suspensions through tubes. I. Single spheres, rods, and discs. Journal of Colloid Science (5), 448-476. 53. Saffman, P. G.; Taylor, G. I., The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences (1242), 312-329. upplementary materials 1 - Area fraction measurements during the imbibition of disk-like and spherical droplets:
Figs. S1 and S2 show how the area fraction of emulsion in field-of-view changes during the imbibition process. This is measured by directly counting the number of the droplets during the imbibition process.
Figure S1. Area fraction of disk-like droplets.
The observed area fraction of disk-like droplets during capillary imbibition for (a) Location 1 and (b) Location 2. The area fraction of each experiment is normalized with the average value observed in that particular experiment. The area fractions stay close to the average values throughout the imbibition. (a) (b)
Figure S2. Area fraction of spherical droplets.
The observed area fraction of disk-like droplets during capillary imbibition for (a) Location 1 and (b) Location 2. The area fraction of each experiment is normalized with the average value observed in that particular experiment. The area fractions stay close to the average values throughout the imbibition. (a) (b)
According to Washburn ’s analysis , during the capillary imbibition of Newtonian fluids, the square of penetration length L(t) is a linear function of time with a slope that depends on the channel and liquid properties. Washburn law is not valid very close to the inlet when the viscous boundary layer is not fully developed, and inertial forces dominate the imbibition process. In the inertial region, the position of meniscus L(t) is a linear function of time . So, the small portion of the experimental data from location one (from the inlet of the channel up to 5 times the width of the channel) is neglected to make sure that the imbibition process is in a viscous-dominated regime. Suspensions are not generally considered Newtonian fluids, but this linear dependence between 𝐿(𝑡) and t is valid, even for concentrated suspension . Motivated by this fact, the Washburn equation's applicability for the case of spherical and disk-like droplets is studied here. Washburn originally studied the imbibition process in circular cross-section channels, but his model can be modified for our system. Since the width of capillaries used in this experiment is 10 times larger than the height, it is reasonable to assume that the velocity gradient in the width direction is negligible with respect to velocity gradient in the height direction and approximation geometry a parallel wall channel. By this assumption, the Naiver-Stokes equation simplifies to: the pressure gradient at time 𝑡 when the liquid has traveled length 𝐿(𝑡) is: Here ∆𝑃 is calculated by the Young-Laplace equation in which 𝜎 is the surface tension, and 𝜃 is the contact angle that is assumed to be constant and equal to zero during the imbibition process. Also, the contribution of meniscus curvature in the y -direction is neglected compared to curvature in the z -direction. Substituting Eq. S2 in Eq. S1 and integrating will result in: 𝜕𝑃𝜕𝑥 = −𝜂 𝜕 𝑣𝜕𝑧 ((S1) 𝜕𝑃𝜕𝑥 = ∆𝑝𝐿(𝑡) = 2𝜎𝐿(𝑡) cos 𝜃 ( 1𝑤 + 1ℎ) ~ 2𝜎𝐿(𝑡)ℎ ((S2) (𝑡) = ℎ𝜎3𝜂 𝑡 = 𝑆𝑡 (S3) Eq. S3 is the modified version of the Washburn equation for the case of parallel wall geometry. Experiments with single-phase liquids confirm the validity of this approximation for the capillary tubes used in this study.
Figure S3 a and b shows the square of penetration length of meniscus normalized by the width of the channel ( 𝐿(𝑡)𝑤 ) as a function of time for different average area fraction of the disk-like emulsion in location 1 and 2. Markers show the experimental measurements and solid lines are the least-squares fit to the data. The dashed lines show the experiment with single phase SDS solutions. The data shows a linear trend for L(t) versus time with 𝑅 > 0.99 for all the experiment. and confirms the validity of Washburn analysis for both locations in the capillary tube, even for the most concentrated emulsion. Also, it confirms that the contact angle is not changing severely during the imbibition. Otherwise, it would appear as a change in the slope during imbibition. Figure S3. Applicability of Wasbhurn to disk-like droplets
The ( 𝐿(𝑡)𝑊 ) versus time for disk-like droplets at (a) location 1 and (b) location 2. The markers show the experimental measurements, the solid-lines are least-squares fit to the data and the dashed-line is the least-squares fit to data from single-phase SDS experimnets. (a) (b) For experiments conducted at location 2, the exact start time of the imbibition is not known. However, it can be estimated by finding the amount of shift it is needed to make the Washburn’s fits pass through the origin. Using this estimation, we can calculate the time meniscus has traveled until it has reached location 2. The same analysis has been conducted for spherical droplets. Results are summarized in Fig. S4:
Figure S4. Applicability of Wasbhurn to disk-like droplets
The ( 𝐿(𝑡)𝑊 ) versus time for spherical droplets at (a) location 1 and (b) location 2. The markers show the experimental measurements, the solid-lines are least-squares fit to the data and the dashed-line is the least-squares fit to data from single-phase SDS experimnets. (a) (b) 𝑳(𝒕) , 𝒍 𝒄 (𝒕) , and 𝒍 𝒔 (𝒕) for disk-like droplets used in the calculation of 𝝂̅ 𝒑𝒑 , 𝝂̅ 𝒕𝒑 and 𝝌 𝒄 : In order to calculate the values of 𝜈 𝑝𝑝 and 𝜈 𝑡𝑝 and 𝜒 𝑐 from capillary imbibition data, we need instantaneous values of 𝐿(𝑡), 𝑙 𝑐 (𝑡), and 𝑙 𝑠 (𝑡). The exact definition of these are depicted in the schematic in supplementary Fig. S5.
Fig. S6. Shows the measured values of
𝐿(𝑡), 𝑙 𝑐 (𝑡), and 𝑙 𝑠 (𝑡) for capillary imbibition of disk-like droplets: Figure S5. Schematic for capillary imbibition of disk-like droplets.
L(t) is the imbibition length at time t, 𝐿 𝑠 (𝑡) is the length of suspension region at time t and 𝐿 𝑐 (𝑡) is the length of clear fluid region. L(t)
Figure S6:
Evolution of 𝐿 ( 𝑡 ) , 𝐿 𝑠 ( 𝑡 ) and 𝐿 𝑐 ( 𝑡 ) versus time for different volume fraction. Top panel is for location 1 and bottom panel is for location 2 - Clausius-Mossotti approximation The expression for Clausius-Mossitti approximation is given as : 𝜈̅ 𝑝𝑝 𝜈̅ = 1 − 12 𝜇̃ 𝜙 𝑎 𝜙 𝑎 (S4) Where 𝜇̃ = 𝑝𝑝 𝜋𝐻 𝑑 (Eq.45, reference 5) and 𝜇 𝑝𝑝 is the polarizability coefficient for a droplet. However, the 𝜇 𝑝𝑝 is not known, but we can estimate an average value from 𝜈̅ 𝑝𝑝 using Eq. 86b from reference 5 𝜈̅ 𝑝𝑝 = 𝐻 𝜈̅ (1 − 12 𝑛̅ 𝑆 𝐻 −3 𝜂 𝜇̅ 𝑝𝑝 ) (S5) Where 𝜇̅ 𝑝𝑝 is the average polarizability coefficient and 𝑛̅ 𝑠 is the number density of the droplets. Substituting the Eq. of S5 in equation S4, results in: 𝜇̃ = 𝐻𝜈̅ − 𝜈̅ 𝑝𝑝 𝑛̅ 𝑠 𝐻𝜈̅ . 1 𝑑 (S6) By substituting the value of 𝜈̅ 𝑝𝑝 and 𝑛̅ 𝑠 from experiment with lowest area fraction, we find 𝜇̃ =3.73 . 𝝌 𝒄 From the plots in Fig. 4a and b, it can be observed that the convective transport coefficient 𝜒 𝑐 shows an increasing trend with respect to the particle density 𝜙̅ 𝑎 . To provide a qualitative explanation we consider a simple model where the particle phase is treated as a porous medium that slides through the channel due to the pressure gradient ∇𝑃 and is restricted by friction from the sidewalls. The fluid phase penetrates through the particle phase with a permeability 𝜇 . The relative velocity of the fluid phase with respect to the particles can be expressed by Darcy’s law as − 𝑈 = −𝜇∇𝑃 (S7) where V and U are the fluid and particle phase velocities. While the balance between the pressure force 𝐹 𝑝 = −ℎ𝑤𝐿∇P and the friction force 𝐹 𝑤𝑓 = 𝑤𝐿𝑓𝑈 responsible for the motion of the porous particle phase gives an expression for the pressure gradient ∇𝑃 = − 𝑓𝑈ℎ (S8) Combining eqn. S7 and S8 we obtain the following relation for relative fluid velocity:
𝑉 − 𝑈 = 𝑓̅𝑈 (S9) where 𝑓 ̅ = 𝑓𝜇/ℎ is the dimensionless friction factor. Assuming a cylindrical shape of droplets due to strong confinement, the velocity of the suspension can be expressed as a weighted sum of the particle phase and fluid phase velocities as 𝑢 = 𝜙̅ 𝑎 U + (1 − 𝜙̅ 𝑎 )𝑉 = 𝑈 + (1 − 𝜙̅ 𝑎 )Δ𝑉 (S10) Using eqn. S9, the above equation can be expressed as a relation between the particle phase velocity u, suspension velocity U, the dimensionless friction factor 𝑓̅ : 𝑢 = 𝑈[1 + (1 − 𝜙̅ 𝑎 )𝑓̅] (S11) Finally, by using eqn. S11 and the definition for convective transport coefficient Eq. 3-3, we arrive at an expression which directly relates the convective transport coefficient 𝜒 𝑐 to the area-fraction of particles 𝜙̅ 𝑎 in the suspension and the dimensionless friction factor 𝑐 = [1 + (1 − 𝜙̅ 𝑎 )𝑓̅] (S12) rom eqn S12 one can observe that as 𝜙̅ 𝑎 → 1 , 𝜒 𝑐 → 1 which emphasizes the impenetrability of the particle phase to fluid flow at very high densities. For intermediate values, one can realize that when area fraction 𝜙̅ 𝑎 increases, the total number of particles in contact with the sidewall also increases leading to an increase in the friction factor f . However, due to the deformability of drops the permeability of the medium falls drastically and at a rate higher than the increase of the friction factor. Additionally, the factor (1 − 𝜙̅ 𝑎 ) also decreases with density. Hence, eqn. S12 decreases strongly with particle density. – Stabilization of the dense-droplet band in stronger confinement
We have conducted experiments in the confinement ratio of 𝑑/ℎ = 0.8 for spherical droplets. We have observed that a dense closely packed droplet band can prevent formation of fingers at this stronger confinement. This is due to the non-zero yield-stress of the ordered emulsion:
Supplementary references
1. Washburn, E. W., The dynamics of capillary flow.
Physical review (3), 273. 2. Quéré, D., Inertial capillarity. EPL (Europhysics Letters) (5), 533. 3. Holloway, W.; Aristoff, J. M.; Stone, H. A., Imbibition of concentrated suspensions in capillaries. Physics of Fluids (8), 081701. 4. Hamraoui, A.; Nylander, T., Analytical approach for the Lucas – Washburn equation.
Journal of colloid and interface science (2), 415-421. 5.
Bławzdziewicz, J.; Wajnryb, E., An analysis of the far -field response to external forcing of a suspension in the Stokes flow in a parallel-wall channel.
Physics of Fluids (9), 093303. Fig. S7 Presence of a closely-packed droplet-band under strong confinements can prevent fingering instability 200 𝛍𝐦𝛍𝐦