The model of a droplet squeezed through a narrow-constricted channel has many applications in pathology, chip/filter/membrane design, drug delivery, etc. Understanding the transient physics of the squeezing process is important in the design and optimization of many micro flow systems. However, available models often ignore the influence of droplet viscosity, and they usually feature low numerical efficiency by solving Navier-Stokes equations. In the present research, we developed a low-dimension analytical model to predict the pressure of squeezing a viscous droplet through a circular constricted channel with acceptable fidelity and low computational cost. Our approach is as follows. We first adapt the Hagen–Poiseuille law to predict the viscosity effect of droplet squeezing. Next, we obtain an analytical expression for the extra pressure caused by only the curvature change obtained. Finally, the general expression of squeezing pressure taking consideration of viscosity and surface tension is expressed. The analytical model we developed is in great agreement with the numerical solutions of the Navier-Stokes equation at a low Reynolds number and low capillary number. These findings have fundamental significance for future applications in engineering and industry.

The model of a single droplet passing through a constricted channel has been widely used in various applications1 such as microfluidics, chemical, biomedical engineering, etc. Specifically, the transient behavior of the droplet squeezing in a constriction plays an important role in cell detection, membrane/filter design, and cell aspiration technologies. Compared to complex particle models (e.g., capsule,2 vesicle,3 compound droplet,4 etc.), a simple droplet model provides quick solutions while staying close to physical reality in many applications, especially for modeling compliant cells and low-viscosity oil droplets.

Numerous research studies have been devoted to revealing the droplet squeezing mechanisms, both experimentally and numerically. For example, recent numerical studies of transient pressure behavior of a droplet squeezing through a constricted channel with various cross sectional geometries5 have identified two types of critical velocities: a co-flow critical velocity6 and a minimum impulse transportation velocity.7 However, insufficient attention has been paid to tackle the problem analytically. Beyond the estimation of critical pressure (by balancing it with surface tension force described by the Young-Laplace equation), we have quite limited knowledge of the transient processes. There are several ways to obtain the pressure profile caused by the interfacial motion during squeezing through a narrow-constricted channel, such as the energy approach,8 and by solving the Navier-Stokes equation analytically (possible for an axially symmetric droplet). Existing analytical models include a droplet squeezed through two parallel plates,9 a bubble squeezed through a short constriction (energy approach),8 and a droplet passing a constriction due to gravity.10 However, further development is inhibited by the complexity of multiphase flow dynamics. Some work involves analyzing and optimizing a single-constricted channel analytically and numerically based on pressure,6 impulse,7 and energy.8 

All the above-mentioned studies are valid in the case when the viscosity of the droplet equals the one of the carrier fluid, i.e., for zero viscosity difference, and there is a limited analytical description of the transient squeezing behavior under the influence of droplet viscosity. To the best of our knowledge, the combined effects of droplet viscosity and droplet curvature change during the squeezing process have not been well investigated. There are only a few experimental11 and numerical5 studies on the dynamics of the viscous droplets squeezing.

However, for many practical applications, the dynamic viscosity of the droplet is much higher than that of the carrier fluid. For example, the viscosity of mammalian cancer cells is approximately 20–40 times of water viscosity.12 The viscosity of oils varies from 10 to 1000 times that of water, and typical biological tissues can reach up to 107 times of water viscosity.13 With further requirements in the processing of highly viscous droplets for biological and chemical applications, such as enhanced oil recovery and tissue aspiration, the transient behavior of viscous droplets becomes even more important. For the systems mentioned above, a more realistic model accounts for both viscosity and surface tension effects.

Recently, researchers have become interested in continuous (constant flow rate boundary conditions) droplet manipulation in a constricted channel, such as droplet generation, cell recognition, drug delivery, etc. For these cases, the transient pressure of a viscous droplet under constant flowrate operation has practical meanings.

In the present research, our motivation is to understand and describe the transient squeezing pressure for a viscous droplet (with viscosity and surface tension) squeezed through a circular narrow-constricted channel with the most efficient analytical treatment. The surface tension term, viscosity term, and background pressure term were derived explicitly. Based on mass conservation, the curvature change has been analytically derived. Furthermore, the Hagen–Poiseuille law was decomposed for a two-phase (droplet and carrier fluid) squeezing model. The combined effect of viscosity and curvature on droplet dynamics has been studied. Finally, the piecewise pressure function was compared with the numerical solution of the Navier-Stokes equation and a good agreement between these two approaches has been demonstrated. These results provided a quick and explicit solution of the viscous droplet squeezing pressure problem that can inspire further research on the complex droplet transportation or constricted channel design.

The physical model of a viscous droplet squeezed through a constricted channel is shown in Fig. 1. The definition of a critical channel length (CCL) and the detailed stages of the squeezing process are given in the  Appendix.

FIG. 1.

Viscous droplet squeezed into a short constricted channel. The droplet is defined by two parameters, surface tension (σ) and viscosity ratio (ηd/ηc).

FIG. 1.

Viscous droplet squeezed into a short constricted channel. The droplet is defined by two parameters, surface tension (σ) and viscosity ratio (ηd/ηc).

Close modal

In the present study, the following assumptions were made:

  • In the curvature change part, the droplet was supposed to be spherical under the action of surface tension or the segment of the sphere under the action of the contact angle.

  • The model was axially symmetric. No off-center cases were considered.14,15

  • The Reynolds number of the carrier fluid was assumed to be low (Re < 1). Because of the typical micrometer scale, the assumption of Stokes flow is always justified.

  • In this study, we used narrow gap assumption with non-dimensionalized gap height, Λ,

Λ=εr1.
(1)

Here, ε is the gap distance between the droplet interface and the constricted channel wall and r is the radius of the constricted channel. The lubrication film thickness is assumed to be much smaller than the constricted channel radius,16 and thus we work under the thin lubrication film assumption. Consequently, we can assume only one phase in the radial direction. The empirical relation shows that the pressure drop induced by the lubrication film can be neglected in a highly viscous droplet at low capillary numbers.17 

  • (e)

    Two key properties of the droplet are surface tension (σ) and viscosity ratio (η*= ηd/ηc) of the droplet (ηd) and the carrier fluid (ηc), respectively. Both surface tension and viscosity require extra driving pressure for a successful squeezing. Other factors such as elasticity, wall-adhesion, viscoelasticity, and bending elasticity of the droplet were neglected.2,18–20

The total flow time of the droplet squeezing is estimated by7 

τ=(VD+VC2V0)/ŪAC.
(2)

Here, VD is the volume of the droplet, VC is the volume of the constricted channel, V0 is the volume of the droplet inside the constricted channel at initial contact, Ū is the average flow velocity in the constricted channel (Ū=1AC0rU2πrdr), and Ac is the area of the constricted channel.

Within the time interval [0, τ], the squeezing pressure is a piecewise function. There are five stages we are interested in (as seen in the  Appendix):

Stage I: 0 < t < T1, the radius of the droplet front cap decreases from R to r. Here, R is the radius of an un-deformed droplet and r is the radius of the constricted channel. The droplet enters the channel with the advancing length varying from h(0) to the critical advancing length hc = r. The advancing length is defined as the droplet penetration depth in the constricted channel, as shown in Fig. 2.

FIG. 2.

Geometry parameters in the droplet squeezing into a constricted channel at initial contact T0.

FIG. 2.

Geometry parameters in the droplet squeezing into a constricted channel at initial contact T0.

Close modal

Due to squeezing at a constant flow rate, T1 reads as

T10=ΔV1ŪAC=1ŪAC12VDV0
(3a)

or

T1=23πr3ŪACV0ŪAC=23rŪV0ŪAC.
(3b)

Here, V0 is the droplet volume inside the constricted channel when the droplet makes an initial contact with the constricted channel.

Stage II: T1 < t < T2 or r < h(t) < L + r, the radius of the front cap stays constant at r, while that of the rear cap continuously decreases. As t = T2, the droplet front cap becomes hemispherical outside the constricted channel so that T2 is given by

T2T1=ΔV2ŪAC=LŪ
(4a)

or

T2=LŪ+23rŪV0ŪAC.
(4b)

Here, the first, the second, and the third terms stand for the flow time inside the constricted channel, the flow time of the droplet hemisphere outside the constricted channel, and the initial volume of the droplet inside the constricted channel. The latter is neglected for simplicity.

Stage III: T2 < t < T3, due to fluid flow, both front and rear curvatures change: the front cap is increasing while the rear cap is shrinking until it reaches a critical exit length (r). Also, there is a moment when the front cap and the rear cap are equal in size.

Stage IV: T3 < t < T4, the flow repeats the inverse process of the droplet entry: the rear cap curvature remains constant, while the front cap outside the exit continuously increases.

Stage V: T4 < t < T5, the radius of the rear cap increases r to its un-deformed radius R, which repeats the inversed process of Stage I.

After T3, the following relations can be obtained from the whole process of the droplet squeezing process ( Appendix):

T3τT2,
(5)
T4τT1,
(6)
T5=τ.
(7)

In general, the source of pressure is either the curvature change or the viscosity length change and the pressure profile can be characterized by five time-pressure points, (Ti, Pi).

Under the assumptions mentioned above, the total squeezing pressure (ΔPtotal) is a function of the Hagen–Poiseuille term inside the constricted channel caused by the carrier fluid (Pvis), the contraction-and-expansion term (ΔPC+E), the surface tension induced pressure term (ΔPσ), and the droplet viscosity induced pressure drop term caused by the viscous droplet (ΔPηd) in the following form:

ΔPtotal=Pvis+ΔPC+E+ΔPσ+ΔPηd.
(8)

The first two terms are constant, collectively called the background pressure drop described in the subsection titled Background pressure drop. In the ΔPσ term, we consider a non-wetting droplet with the contact angle of 180°.

The Hagen term and the contraction-and-expansion term are collectively called the background pressure. The Hagen law viscosity-induced pressure drop is given by

Pvis=8η0Lr2Ū.
(9)

Here, η0 is the dynamic viscosity of the carrier fluid, r is the radius of the constricted channel, L is the length of the constricted channel, and Ū is the average carrier fluid velocity in the constricted channel.

The minor loss, corresponding to the second-order velocity effect, is

PCE=KE+KCρŪ22,
(10)

where KE + KC = 1.5, for a sudden constriction from a large to small area.

For a semi-static operation condition, the curvature-induced Young-Laplace pressure drop can be expressed as a function of two curvatures, 1/RR (right-side droplet cap) and 1/RL (left-side droplet cap),

Δpσ=2σ(1/RR1/RL).
(11)

A detailed procedure to obtain Laplace pressure is described for Stage I.

As seen in Fig. 2, h(t) reaches r at T1. Thus, κ(T1) = 1/r, and the system operates under a critical flow condition with the highest interface resistance.

The curvature of the droplet spherical cap is a function of h(t) and r,

κ(t)=2h(t)h2(t)+r2.
(12)

Here, h(t) is the height of the spherical cap at flow time t. In Stage I, h(t) is the advancing length inside the constricted channel. At initial contact (the droplet contacts with the channel wall), the advancing length inside the constricted channel is

h(0)=RR2r2.
(13)

Both the left and right caps are assumed to be segments of a sphere. The volume change of the front cap can be obtained by integrating the continuity equation dV=ŪAcdt or

VR(t)=0tŪAcdt+V0.
(14)

Here, V0=π6h0r2+h2(0), with h(0), Ū, and Ac being the initial advancing length, the average flow velocity in the constricted channel (Ū=1AC0rU2πrdr), and the area of the constricted channel (Ac = πr2), respectively, and VD and VC stand for the volume of the droplet and the volume of the constricted channel (VC = AcL).

Thus, for distances smaller than the critical advancing length (h(t) < r), the flow is purely curvature-based. The volume of the droplet inside the constricted channel can be expressed as Vt=π6h(t)(3r2+h2(t)). Combining the above equations, the relation between the volume of the droplet in the constricted channel and the advancing length at time t reads as

π6h3(t)+π2r2h(t)Vt=0
(15a)

or

h3t+3r2h(t)6πVt=0.
(15b)

If written in a general form, this equation reads as x3 + px + q = 0, with the analytical solution obtained using the Cardano method x=q2+q22+p333+q2q22+p333, where p = 3r2 and q=6πV(t). Then, the dynamic advancing length inside the constricted channel can be derived from

h(t)=3πVt+3πVt2+r63+3πVt3πVt2+r63,
(16)

where Vt=ŪAct+V0=ŪAct+π6h0r2+h02. For different stages, V0 and V(t) vary.

For other stages, a similar approach can be used to get the curvature through the height of the cap and the squeezing time.

In the case of surface tension dominating, we assume the viscosity ratio to be η* = 1. For other cases with η* ≠ 1, a two-phase system is formed along the channel direction, which consists of the droplet with the viscosity ηd and the carrier fluid with viscosity η0. The viscosity-induced pressure drops for different stages are presented in the following paragraphs.

Stage I: There is no effective contact length at the channel wall. Thus, the viscosity-induced pressure drop is 0.

Stage II: The effective contact length increases linearly with flow time L1=ŪΔt. Here, Δt is defined as Δt = t − T1. Thus, the viscosity-induced pressure drop for the left-side (ΔPL) of the viscous droplet is given by

ΔPL=8r2ηdŪΔtŪ.
(17a)

The viscosity-induced pressure drop for the right side (ΔPR) in the carrier fluid reads as

ΔPR=8r2η0LŪΔtŪ,
(17b)

where Δt is the flow time, effective within the time domain, Δt > T1. T1 was defined in the time analysis part, in Eq. (3b). Applying the superposition principle, we arrive at

ΔPη=ΔPL+ΔPR=8r2ηdŪΔtŪ+8r2η0LŪΔtŪ.
(18a)

Here subscripts L and R stand for the fluid left side of interface and right side of the droplet front cap interface, respectively. Rearranging the terms with respect to the flow time, we get

ΔPη=ΔPη(t)+ΔPη(const)=8r2(ηdη0)Ū2Δt+8r2η0LŪ.
(18b)

The above equation is the time-dependent pressure drop for the viscous droplet entering the constricted channel. As can be seen from Eq. (18b), the transient pressure of the viscous droplet increases linearly with flow time. The slope can be determined from k=8r2(ηdη0)Ū2, which expresses the proportionality of the slope to two variables, the viscosity difference ηdη0, and the velocity squared, Ū2. Besides, the pressure also needs to overcome the constant term ΔPμ(const)=8r2η0LŪ, that is, the background pressure with no droplet formed, according to the single-phase Hagen-Poiseuille equation.

The chamber stands for the unconstricted channel at the upstream and downstream of the constricted channel. The pressure drop in the chamber is much smaller than the one in the constricted channel, Pchamber1 = Pchamber2 ≈ 0. The slope of the P-L curve (Fig. 3) changes with the viscosity of the droplet. With the viscous droplet advancing, the pressure increases from Pin to Pin′. At the inlet, the pressure drop is proportional to the contact length between the viscous droplet and the channel wall.

FIG. 3.

The pressure drop for a pure viscous droplet squeezed into the constricted channel (P-L curve). The shaded area is the constricted channel. L1 is the effective contact length. The slopes of the pressure droplet are different due to the viscosity difference.

FIG. 3.

The pressure drop for a pure viscous droplet squeezed into the constricted channel (P-L curve). The shaded area is the constricted channel. L1 is the effective contact length. The slopes of the pressure droplet are different due to the viscosity difference.

Close modal

Stage III: The droplet filled the channel, and the contact length stays constant (L). The viscosity-induced pressure drop is a constant (8ηdLr2Ū).

Stage IV: The droplet begins to pull out of the constricted channel. Consequently, L1 is occupied by the carrier fluid, while L2 is occupied by the viscous droplet. Thus,

ΔPη=ΔPL+ΔPR=8r2η0ŪΔtŪ+8r2ηdLŪΔtŪ
(19a)

or

ΔPη=ΔPηt+ΔPηconst=8r2(ηdη0)Ū2Δt+8r2ηdLŪ.
(19b)

Here, Δt = t − T3.

Stage V repeats Stage I, and there is no droplet viscosity-induced pressure drop.

According to the squeezing pressure expression derived above, for a pure viscous droplet model, there are three pressure regimes, as seen in Fig. 4: (1) a square pressure signal for an extremely short channel, (2) a trapezoidal pressure signal for a short-constricted channel, and (3) a triangle pressure signal in a channel at CCL.

FIG. 4.

Three types of pressure for a pure viscous droplet for different channel lengths. (a) Entry length is extremely short (square), (b) a normal short constricted channel (trapezoid), and (c) a constricted channel (triangle) with the length of CCL, as defined in the  Appendix.

FIG. 4.

Three types of pressure for a pure viscous droplet for different channel lengths. (a) Entry length is extremely short (square), (b) a normal short constricted channel (trapezoid), and (c) a constricted channel (triangle) with the length of CCL, as defined in the  Appendix.

Close modal

From the above analysis, we found that the total pressure drop is determined by the background pressure (Pb) as well as the following three variables: (1) the curvature of the droplet inside the constricted channel, (2) the curvature of the droplet outside the constricted channel, and (3) the contact length within the constricted channel. Figure 5 show the pressure profile and the droplet deformation of a two-component droplet in (a) a pure viscous droplet, (b) a pure surface tension droplet, and (c) the combination of two components.

FIG. 5.

Transit pressure model for a droplet with (a) pure viscosity, (b) pure surface tension, and (c) a combination of viscosity and surface tension, entering a circular constricted channel at a constant velocity.

FIG. 5.

Transit pressure model for a droplet with (a) pure viscosity, (b) pure surface tension, and (c) a combination of viscosity and surface tension, entering a circular constricted channel at a constant velocity.

Close modal

The Volume of Fluid (VOF) algorithm is employed to capture the interface between the droplet and the carrier fluid.21 This algorithm is computationally friendly with low computer storage requirements. It has been implemented in many solvers such as Fluent, OpenForm, Polyflow, and Star-CCM. In this study, we used a commercial software, ANSYS Fluent, which is the most popular fluid flow solver providing robust and quick solutions.

The interface is tracked by solving the transport equation for the primary-phase volume fraction function. Assuming no mass transfer between the two phases, the continuity equation can be written as

αt+(αU)=0.
(20)

Here U is the velocity vector. The substantial derivative of the primary-phase volume fraction is zero so that the interface-tracking is achieved through the velocity field at the interface.

The momentum equation is solved for the whole computational domain, and the velocity filed is shared between the two phases,

t(ρU)+(ρUU)=p+μU+UT+ρg+FS,
(21)

where g is the gravitational acceleration vector and FS is the surface tension force source. ρ and μ are the volume-fraction-averaged density and viscosity, respectively.

The volume of fluid (VOF) method has been widely used in the simulation of a microfluidics system.22 The VOF algorithm assumes a single-fluid problem formulation with the volume-fraction-averaged material properties at the finite-thickness interface. It employs a volume fraction function (α) of the primary phase (the carrier fluid), which varies between 0 (droplet) and 1 (main phase). An explicit time stepping scheme is used with the Courant number within 0.5-5 to ensure the stability and convergence. The cell-water interface is reconstructed using the Geo-reconstruct scheme available in Fluent. Inflation layers are used to achieve an accurate resolution of the boundary layer. The mesh is inflated on the wall of the narrow-constricted channel.

Another challenge of numerical simulation is the discontinuity in material properties (such as viscosity) at the interface. Those properties are derived from the volume fraction average of the droplet and the carrier fluid. For a two-phase flow system, the volume-fraction-averaged viscosity can be calculated as follows:

η=αηc+(1α)ηd.
(22)

Here ηc and ηd stand for the viscosity of the carrier fluid and the droplet, respectively. Neglecting the influence of gravity, we assume the densities of two phases to be the same, i.e., ρ = ρc = ρd. These properties are then used to get the velocity field from a single momentum equation in the domain. There is no explicit limit of the viscosity value in the Ansys Fluent solver. However, from the Carreau model (a non-Newtonian fluid), it can be seen that the viscosity setting can reach up to 14.92 Pa s as an upper limit and 6.433 mPa s as a lower limit. Similarly, we use the viscosity ratio as a modeling limit with 14.92 Pa s as a fluid upper limit.23 

The Laplace pressure can be expressed by the curvature change, κ=n,

Δp=κσ=σn,
(23)

where σ is the interfacial tension coefficient between the droplet and the carrier fluid and κ is the mean curvature of the interface. The continuum surface force (CSF) model, proposed by Brackbill et al.,24 was used to simulate the interfacial-tension-induced pressure jump across the interface. In this method, the interfacial surface force is transformed into a volume force using the divergence theorem. In the limit of zero interface thickness, the interfacial tension term in the volume force reads as

FS=σ2ρκα(ρ1+ρ2)=σκα.
(24)

This interfacial tension force is the source term in the momentum equation Eq. (21).

In a narrow confinement, the droplet interacts with the constricted channel wall. The wall adhesion boundary condition in the rigid constricted channel can be obtained by imposing a dynamic boundary condition and considering the effect of the contact angle with the unit normal vector

ni=nwcosθd+ntsin(θd),
(25)

where nw is the unit vector normal to the wall, nt is a vector tangent to the wall and normal to the contact line, and θd is the dynamic contact angle between the fluid and the wall. θd depends on the fluid properties, smoothness, and geometry of the wall, as well as the flow velocity in the constricted channel. In this study, the dynamic contact angle was set at 180° for the calculation of the wall adhesion.

Due to axial-symmetry, the boundary condition at the center of the constricted channel is ur = 0, so a fully developed inlet velocity profile is parabolic. The wall is considered impermeable and rigid. A non-slip boundary condition is assumed at the wall, ui,wall = 0. The pressure at the outlet is set as a reference (zero-gauge pressure), Pout = 0.

The droplet front cap deforms to achieve more acute curvature and adjust an input pressure. The droplet may be balanced by the static friction at the tip. Beyond the critical point, the droplet interacts with the channel wall and the curvature recovers to a hemisphere. This process was described in previous studies.6 

Stage I was implemented using Eqs. (11), (12), and (16). For other stages, the curvature expression and height of the spherical cap are changed accordingly.

In Stage VI and Stage V, the droplet recovers. The surface tension induced pressure drop recovers in the following mathematical form for T3 < t < T5:

Pσt+Pστt=28r2ηdLŪ.
(26)

Thus, there is a moment of time (Tm = τ/2), after which the surface tension induced pressure becomes symmetric about flow time (Tm) and pressure (ΔPμ=8r2ηdLŪ).

Geometrically speaking, the viscosity-induced squeezing pressure is symmetric at t = Tm in the following form:

Pvist+2σ1r1R=Pvisτt2σ(1r1R).
(27)

The numerical model has been validated in our previous work.5–7 More validations are not discussed in the present work.

To reduce the number of variables, we constructed a group of dimensionless parameters combining the variables given in Table I.

TABLE I.

Variables and corresponding dimensions.

VariableDescriptionsDimensions
Squeezing pressure ML−1T−2 
Channel length 
Radius of an undeformed droplet 
Radius of the constricted channel 
Σ Surface tension MT−2 
Fluid density ML−3 
Flow time 
Ū Average velocity in the constricted channel LT−1 
η0 Dynamic viscosity of carrier fluid ML−1T−1 
ηd Dynamic viscosity of the droplet ML−1T−1 
VariableDescriptionsDimensions
Squeezing pressure ML−1T−2 
Channel length 
Radius of an undeformed droplet 
Radius of the constricted channel 
Σ Surface tension MT−2 
Fluid density ML−3 
Flow time 
Ū Average velocity in the constricted channel LT−1 
η0 Dynamic viscosity of carrier fluid ML−1T−1 
ηd Dynamic viscosity of the droplet ML−1T−1 

As outlined in the literature, the squeezing pressure is a function of nine independent variables

Pt=fL,R,r,σ,ρ,Ū,ηd,η0,t.
(28)

The density difference is neglected since we do not consider gravity in this study. From Table I, we can also reveal that all ten variables are formed via the combination of units of mass (M), length (L), and time (T). We use the Buckingham Pi theorem to perform the non-dimensionalization. We construct Pi groups for the squeezing pressure first since it is the dependent variable

Π1=Praσbρc,
(29)
{M0L0T0}=(MLT2)(L)a(MT2)b(ML3)c.
(30)

Equating the powers of mass, length, and time, we get

1+b+c=0,a13c=0b+1=0.,
(31)

With a = 1, b = −1, and c = 0, we arrive at

Π1=Ptσ/r.
(32a)

However, to get a meaningful non-dimensional parameter, we divide it by the critical pressure (derived from the Young- Laplace equation)

Π1=P*=Pt2σ(1r1R).
(32b)

Next, two length parameters were non-dimensionalized by the radius of the constricted channel (r),

Π2=R*=Rr,
(33)
Π3=L*=Lcr.
(34)

The viscosity of the droplet was non-dimensionalized by the viscosity of the carrier fluid,

Π4=η*=ηdη0.
(35)

The transient time was non-dimensionalized with the transient time,

Π5=t*=tτ.
(36)

The surface tension was non-dimensionalized as

Π6=σρarbŪc,
(37a)

where a = 0, b = −1, and c = 1. Introducing the capillary number (Ca), we get

Π6=Ca=μ0Ūσ.
(37b)

Here, we take the ratio of Hagen-Poiseuille’s term (8ηdLr2Ū) and Young-Laplace (2σr) one and assume that the droplet radius is much bigger than that of the constricted channel, 8ηdLŪr2/2σr=4LηdŪrσ=4LrCad. This ratio tells us that the dominating force (capillary or viscous) is determined by the non-dimensional channel length and the droplet-based capillary number.

The viscosity term was non-dimensionalized using the Reynolds number,

Π7=Re=2ρŪrμ0.
(38)

Finally, after the non-dimensionalization of the variables, we arrive at

Pt2σ(1r1R)=f{Lr,Rr,ηdη0,tτ,Ca,Re}.
(39)

In this study, we fixed the constricted channel radius (r = 2.5 µm), droplet radius (R/r = 3.2), and density of both fluids (ρ = 103 kg/m3).

Figure 6 is the comparison of the analytical model and numerical results for the non-viscous droplet with different surface tension coefficients. It can be seen that the analytical model agrees well with the numerical data (with some fluctuations due to numerical errors).

FIG. 6.

Comparisons of numerical and analytical results for two droplets with a surface tension of 50 and 10 mN/m, respectively, at Re = 0.18 for (a) half time in an x-log coordinate and (b) the whole process.

FIG. 6.

Comparisons of numerical and analytical results for two droplets with a surface tension of 50 and 10 mN/m, respectively, at Re = 0.18 for (a) half time in an x-log coordinate and (b) the whole process.

Close modal

Figure 7 shows the comparisons between analytical and numerical models for a pure viscous droplet. At high viscosity, the droplet might breakup at the exit. The rectangle signal for the pressure is not provided since it is the case of an extremely short channel with high viscosity, which is not stable numerically.

FIG. 7.

Comparision of analytical and numerical results for a high viscous droplet. (a) Pressure profile in a short channel. (b) A constricted channel with the length of CCL.

FIG. 7.

Comparision of analytical and numerical results for a high viscous droplet. (a) Pressure profile in a short channel. (b) A constricted channel with the length of CCL.

Close modal

Figure 8 is a comparison of analytical and numerical results for a droplet taking consideration of both surface tension and viscosity ratio (η* = 10). At a low Reynolds number (<1), the analytical solutions match well with the Navier-Stokes solutions.

FIG. 8.

Comparision of the analytical model with numerical results for a droplet with a viscosity ratio of 10.

FIG. 8.

Comparision of the analytical model with numerical results for a droplet with a viscosity ratio of 10.

Close modal

As the Reynolds number increases, the curvature of the front cap increases as well. Before co-flow happens, the curvature increase will increase the squeezing pressure. The numerical results obtained by solving the Navier-Stokes equations are higher than the corresponding analytical predictions with a spherical assumption. Besides, the slop of the pressure drop in numerical results is smaller than the pure viscosity model prediction because the pure viscosity model neglected the curvature change. However, with the decrease of the surface tension (a lower capillary number), the prediction becomes more accurate (Fig. 9).

FIG. 9.

Pressure predictions for different viscosities at Re = 1.8.

FIG. 9.

Pressure predictions for different viscosities at Re = 1.8.

Close modal

The model we developed in the present research demonstrated excellent accuracy in a short-constricted channel at a low Reynolds number and low capillary number. However, this model has not taken consideration of the dynamic effect and the complicated droplet deformation at the channel exit. Therefore, the discrepancy was shown in the higher Reynolds number, high capillary number, or long-constricted channel, which is beyond the assumption of the research. Figure 10 shows three sources causing prediction discrepancy: (a) co-flow at the inlet, (b) dynamic effect at the outlet, and (c) un-uniform velocity field due to the interface.

FIG. 10.

Sources of inaccuracy. Reasons for the discrepancy between numerical and the analytical solutions: (a) co-flow at a high Reynolds number (>1), (b) spherical cap assumption failure due to the high Reynolds number and high viscosity ratio, and (c) the velocity field inside the constricted channel is not uniform due to the existence of the interface and viscosity difference.

FIG. 10.

Sources of inaccuracy. Reasons for the discrepancy between numerical and the analytical solutions: (a) co-flow at a high Reynolds number (>1), (b) spherical cap assumption failure due to the high Reynolds number and high viscosity ratio, and (c) the velocity field inside the constricted channel is not uniform due to the existence of the interface and viscosity difference.

Close modal

In addition, the breakup of a viscous droplet, which was found to be related to the degree of confinement, viscosity ratio, and Ca number,25,26 has not been taken into consideration. Due to the complexity of the instability problem, the influence of the droplet breakup on the squeezing pressure has not been reported and considered. Briefly speaking, the expansion of the droplet at the outlet as described above, leading to a sudden decrease of the local pressure, caused the breakup of the droplet.27 

In the present research, we analytically expressed the squeezing pressure of a viscous droplet squeezing through a circular constricted channel. Two properties, surface tension and viscosity ratio, were taken into consideration. The transient squeezing pressure was found to be a piece-wise function determined by two variables: curvature and advancing length. By comparison with numerical results, the analytical model demonstrated high accuracy under five assumptions. Furthermore, we noticed the limitation of the present model to predict high-velocity cases. The model can be improved and extended by (1) evaluating the analytical model quantitatively by more cases through larger samples in different applications, (2) considering flow instability such as the droplet breakup in the study of printing or droplet generation, and (3) adding complex constitutive equations, such as viscoelastic properties, to simulate applications such as cancer metastasis or drug delivery.

Zhifeng Zhang thanks the support of Harvey & Geraldine Brush Fellowship and Max & Joan Schlienger Scholarship awarded by the College of Engineering, the Pennsylvania State University. The project is partially supported by NSF Grant (No. 1711798) awarded to Jie Xu, at the University of Illinois at Chicago.

A long-constricted channel is defined as a channel that can squeeze the whole droplet inside the channel. Then the critical channel length (CCL) is the one, where the droplet can be fully squeezed (Fig. 11). When the constricted channel is shorter than CCL, the curvature changes at both sides of it (Fig. 12). If the channel is longer than CCL, there is only one-side curvature change and the curvature of the other side of the cap remains constant.

FIG. 11.

(a) Definition of the critical length and relationship between critical length and droplet radius and (b) channel radius relationship for (I) short channel; (II) channel with critical length; and (III) channel longer than the critical length.

FIG. 11.

(a) Definition of the critical length and relationship between critical length and droplet radius and (b) channel radius relationship for (I) short channel; (II) channel with critical length; and (III) channel longer than the critical length.

Close modal
FIG. 12.

Squeezing stage difference in (a) a short-constricted channel and (b) a long-constricted channel.

FIG. 12.

Squeezing stage difference in (a) a short-constricted channel and (b) a long-constricted channel.

Close modal

CCL is determined by the radius of the droplet and the radius of the constricted channel, VD = πr2Lc + Vd. Here, Vd=43πr3. The relationship between L/r and R/r in the non-dimensional form is

CCLr=43(Rr)31.
(A1)
1.
Z.
Zhang
,
J.
Xu
, and
C.
Drapaca
, “
Particle squeezing in narrow confinements
,”
Microfluid. Nanofluid.
22
,
120
(
2018
).
2.
Z. Y.
Luo
,
L.
He
, and
B. F.
Bai
, “
Deformation of spherical compound capsules in simple shear flow
,”
J. Fluid Mech.
775
,
77
(
2015
).
3.
E.
Benet
and
F.
Vernerey
, “
Mechanics and stability of vesicles and droplets in confined spaces
,”
Phys. Rev. E
94
,
062613
(
2016
).
4.
Z.
Zhang
,
J.
Xu
, and
X.
Chen
,
Compound Droplet Modelling of Circulating Tumor Cell Microfiltration
(
American Society of Mechanical Engineers
,
2015
).
5.
Z.
Zhang
,
J.
Xu
,
B.
Hong
, and
X.
Chen
, “
The effects of 3D channel geometry on CTC passing pressure–towards deformability-based cancer cell separation
,”
Lab Chip
14
,
2576
(
2014
).
6.
Z.
Zhang
,
X.
Chen
, and
J.
Xu
, “
Entry effects of droplet in a micro confinement: Implications for deformation-based circulating tumor cell microfiltration
,”
Biomicrofluidics
9
,
024108
(
2015
).
7.
Z.
Zhang
,
C.
Drapaca
,
X.
Chen
, and
J.
Xu
, “
Droplet squeezing through a narrow constriction: Minimum impulse and critical velocity
,”
Phys. Fluids
29
,
072102
(
2017
).
8.
M. J.
Jensen
,
G.
Goranović
, and
H.
Bruus
, “
The clogging pressure of bubbles in hydrophilic microchannel contractions
,”
J. Micromech. Microeng.
14
,
876
(
2004
).
9.
D.
Barthès-Biesel
,
Microhydrodynamics and Complex Fluids
(
CRC Press
,
2012
).
10.
A.
Marmur
, “
Penetration of a small drop into a capillary
,”
J. Colloid Interface Sci.
122
,
209
(
1988
).
11.
S.
Ma
,
J. M.
Sherwood
,
W. T. S.
Huck
, and
S.
Balabani
, “
On the flow topology inside droplets moving in rectangular microchannels
,”
Lab Chip
14
,
3611
(
2014
).
12.
T.
Kalwarczyk
,
N.
Ziebacz
,
A.
Bielejewska
,
E.
Zaboklicka
,
K.
Koynov
,
J.
Szymanski
,
A.
Wilk
,
A.
Patkowski
,
J.
Gapinski
,
H.-J.
Butt
, and
R.
Holyst
, “
Comparative analysis of viscosity of complex liquids and cytoplasm of mammalian cells at the nanoscale
,”
Nano Lett.
11
,
2157
(
2011
).
13.
K.
Guevorkian
,
M.-J.
Colbert
,
M.
Durth
,
S.
Dufour
, and
F.
Brochard-Wyart
, “
Aspiration of biological viscoelastic drops
,”
Phys. Rev. Lett.
104
,
218101
(
2010
).
14.
Y.
Ma
,
M.
Zheng
,
M. G.
Bah
, and
J.
Wang
, “
Effects of obstacle lengths on the asymmetric breakup of a droplet in a straight microchannel
,”
Chem. Eng. Sci.
179
,
104
(
2018
).
15.
Z.
Luo
and
B.
Bai
, “
Off-center motion of a trapped elastic capsule in a microfluidic channel with a narrow constriction
,”
Soft Matter
13
,
8281
(
2017
).
16.
F.
Bretherton
, “
The motion of long bubbles in tubes
,”
J. Fluid Mech.
10
,
166
(
1961
).
17.
Q.
Li
and
P.
Angeli
, “
Experimental and numerical hydrodynamic studies of ionic liquid-aqueous plug flow in small channels
,”
Chem. Eng. J.
328
,
717
(
2017
).
18.
M.
Nooranidoost
,
D.
Izbassarov
, and
M.
Muradoglu
, “
Droplet formation in a flow focusing configuration: Effects of viscoelasticity
,”
Phys. Fluids
28
,
123102
(
2016
).
19.
S. G.
Sontti
and
A.
Atta
, “
CFD analysis of microfluidic droplet formation in non-Newtonian liquid
,”
Chem. Eng. J.
330
,
245
(
2017
).
20.
Y.
Wang
,
M.
Do-Quang
, and
G.
Amberg
, “
Viscoelastic droplet dynamics in a Y-shaped capillary channel
,”
Phys. Fluids
28
,
033103
(
2016
).
21.
T.
Darvishzadeh
and
N. V.
Priezjev
, “
Effects of crossflow velocity and transmembrane pressure on microfiltration of oil-in-water emulsions
,”
J. Membr. Sci.
423-424
,
468
(
2012
).
22.
M.
Nekouei
and
S.
Vanapalli
, “
Volume-of-fluid simulations in microfluidic T-junction devices: Influence of viscosity ratio on droplet size
,”
Phys. Fluids
29
,
032007
(
2017
).
23.
C.
Rodriguez-Rivero
,
E. M. M.
del Valle
, and
M. A.
Galán
, “
CFD study of capillary jets under superimposed destabilizing conditions for microdroplet formation
,”
Eng. Appl. Comput. Fluid Mech.
9
,
419
(
2015
).
24.
J.
Brackbill
,
D. B.
Kothe
, and
C.
Zemach
, “
A continuum method for modeling surface tension
,”
J. Comput. Phys.
100
,
335
(
1992
).
25.
N.
Barai
and
N.
Mandal
, “
Breakup modes of fluid drops in confined shear flows
,”
Phys. Fluids
28
,
073302
(
2016
).
26.
S.
Guido
and
V.
Preziosi
, “
Droplet deformation under confined Poiseuille flow
,”
Adv. Colloid Interface Sci.
161
,
89
(
2010
).
27.
L.
Shui
,
F.
Mugele
,
A.
van den Berg
, and
J.
Eijkel
, “
Geometry-controlled droplet generation in head-on microfluidic devices
,”
Appl. Phys. Lett.
93
,
153113
(
2008
).