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.

## INTRODUCTION

The model of a single droplet passing through a constricted channel has been widely used in various applications^{1} 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 geometries^{5} have identified two types of critical velocities: a co-flow critical velocity^{6} 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 experimental^{11} and numerical^{5} 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 10^{7} 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.

## PHYSICAL MODEL AND ASSUMPTIONS

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.

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, Λ,

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}

## ANALYTICAL METHOD

### Flow time analysis

The total flow time of the droplet squeezing is estimated by^{7}

Here, *V*_{D} is the volume of the droplet, *V*_{C} is the volume of the constricted channel, *V*_{0} is the volume of the droplet inside the constricted channel at initial contact, $\u016a$ is the average flow velocity in the constricted channel ($\u016a=1AC\u222b0rU2\pi rdr$), and *A*_{c} 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 < *T*_{1}, 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 h_{c} = r. The advancing length is defined as the droplet penetration depth in the constricted channel, as shown in Fig. 2.

Due to squeezing at a constant flow rate, *T*_{1} reads as

or

Here, V_{0} is the droplet volume inside the constricted channel when the droplet makes an initial contact with the constricted channel.

Stage II: *T*_{1} < t < *T*_{2} 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 = *T*_{2}, the droplet front cap becomes hemispherical outside the constricted channel so that *T*_{2} is given by

or

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: T_{2} < t < T_{3}, 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: T_{3} < t < T_{4}, 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: T_{4} < t < T_{5}, the radius of the rear cap increases r to its un-deformed radius R, which repeats the inversed process of Stage I.

After T_{3}, the following relations can be obtained from the whole process of the droplet squeezing process ( Appendix):

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, (*T*_{i}, *P*_{i}).

## PRESSURE MODEL

### Total pressure drop

Under the assumptions mentioned above, the total squeezing pressure ($\Delta Ptotal$) is a function of the Hagen–Poiseuille term inside the constricted channel caused by the carrier fluid (*P*_{vis}), the contraction-and-expansion term ($\Delta PC+E$), the surface tension induced pressure term ($\Delta P\sigma $), and the droplet viscosity induced pressure drop term caused by the viscous droplet ($\Delta P\eta \u2212d$) in the following form:

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

### Background pressure drop

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

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 $\u016a$ is the average carrier fluid velocity in the constricted channel.

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

where *K*_{E} + *K*_{C} = 1.5, for a sudden constriction from a large to small area.

### Curvature-induced pressure drop

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

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

As seen in Fig. 2, *h*(t) reaches *r* at T_{1}. Thus, *κ*(*T*_{1}) = 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*,

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

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=\u016aAcdt$ or

Here, $V0=\pi 6h0r2+h2(0)$, with *h*(0), $\u016a$, and *A*_{c} being the initial advancing length, the average flow velocity in the constricted channel ($\u016a=1AC\u222b0rU2\pi rdr$), and the area of the constricted channel (*A*_{c} = *πr*^{2}), respectively, and *V*_{D} and *V*_{C} stand for the volume of the droplet and the volume of the constricted channel (*V*_{C} = *A*_{c}*L*).

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 V$t=\pi 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

or

If written in a general form, this equation reads as *x*^{3} + *px* + *q* = 0, with the analytical solution obtained using the Cardano method $x=\u2212q2+q22+p333+\u2212q2\u2212q22+p333$, where *p* = 3*r*^{2} and $q=\u22126\pi V(t)$. Then, the dynamic advancing length inside the constricted channel can be derived from

where $Vt=\u016aAct+V0=\u016aAct+\pi 6h0r2+h02$. For different stages, *V*_{0} 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.

### Viscosity-induced pressure drop

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=\u016a\Delta t$. Here, Δt is defined as Δt = t − *T*_{1}. Thus, the viscosity-induced pressure drop for the left-side ($\Delta PL$) of the viscous droplet is given by

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

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

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

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(\eta d\u2212\eta 0)\u016a2$, which expresses the proportionality of the slope to two variables, the viscosity difference *η*_{d} − *η*_{0}, and the velocity squared, $\u016a2$. Besides, the pressure also needs to overcome the constant term $\Delta P\mu (const)=8r2\eta 0L\u016a$, 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, *P*_{chamber1} = *P*_{chamber2} ≈ 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 *P*_{in} to *P*_{in}′. At the inlet, the pressure drop is proportional to the contact length between the viscous droplet and the channel wall.

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

Stage IV: The droplet begins to pull out of the constricted channel. Consequently, L_{1} is occupied by the carrier fluid, while L_{2} is occupied by the viscous droplet. Thus,

or

Here, Δt = t − *T*_{3}.

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.

### Pressure combination and evolution

From the above analysis, we found that the total pressure drop is determined by the background pressure (*P*_{b}) 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.

## NUMERICAL METHOD

### Volume of Fluid (VOF) method

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

Here $U\u21c0$ 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,

where $g\u21c0$ is the gravitational acceleration vector and $F\u21c0S$ 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.

### Material modeling

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:

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}

### Continuum surface force (CSF) method

The Laplace pressure can be expressed by the curvature change, $\kappa =\u2207\u22c5n$,

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

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

### Boundary conditions

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

where $n\u21c0w$ is the unit vector normal to the wall, $n\u21c0t$ 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 *u*_{r} = 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, *u*_{i,wall} = 0. The pressure at the outlet is set as a reference (zero-gauge pressure), *P*_{out} = 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}

## RESULTS

### Model implementation

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 *T*_{3} < *t* < *T*_{5}:

Thus, there is a moment of time (*T*_{m} = τ/2), after which the surface tension induced pressure becomes symmetric about flow time (*T*_{m}) and pressure ($\Delta P\mu =8r2\eta dL\u016a$).

Geometrically speaking, the viscosity-induced squeezing pressure is symmetric at *t* = *T*_{m} in the following form:

### Non-dimensionalization

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.

Variable . | Descriptions . | Dimensions . |
---|---|---|

P | Squeezing pressure | ML^{−1}T^{−2} |

L | Channel length | L |

R | Radius of an undeformed droplet | L |

R | Radius of the constricted channel | L |

Σ | Surface tension | MT^{−2} |

P | Fluid density | ML^{−3} |

T | Flow time | T |

$\u016a$ | Average velocity in the constricted channel | LT^{−1} |

η_{0} | Dynamic viscosity of carrier fluid | ML^{−1}T^{−1} |

η_{d} | Dynamic viscosity of the droplet | ML^{−1}T^{−1} |

Variable . | Descriptions . | Dimensions . |
---|---|---|

P | Squeezing pressure | ML^{−1}T^{−2} |

L | Channel length | L |

R | Radius of an undeformed droplet | L |

R | Radius of the constricted channel | L |

Σ | Surface tension | MT^{−2} |

P | Fluid density | ML^{−3} |

T | Flow time | T |

$\u016a$ | Average velocity in the constricted channel | LT^{−1} |

η_{0} | Dynamic viscosity of carrier fluid | ML^{−1}T^{−1} |

η_{d} | Dynamic viscosity of the droplet | ML^{−1}T^{−1} |

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

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

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

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

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

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

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

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

The surface tension was non-dimensionalized as

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

Here, we take the ratio of Hagen-Poiseuille’s term ($8\eta dLr2\u016a$) and Young-Laplace $(2\sigma r)$ one and assume that the droplet radius is much bigger than that of the constricted channel, $8\eta dL\u016ar2/2\sigma r=4L\eta d\u016ar\sigma =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,

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

In this study, we fixed the constricted channel radius (r = 2.5 *µ*m), droplet radius (R/r = 3.2), and density of both fluids (*ρ* = 10^{3} kg/m^{3}).

### Numerical validation

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).

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.

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.

### Reynolds number

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).

### Model limitation and error analysis

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.

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}

## CONCLUSIONS

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.

## ACKNOWLEDGMENTS

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.

### APPENDIX: CRITICAL CHANNEL LENGTH (CCL)

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.

CCL is determined by the radius of the droplet and the radius of the constricted channel, *V*_{D} = *πr*^{2}*L*_{c} + *V*_{d}. Here, $Vd=43\pi r3$. The relationship between L/r and R/r in the non-dimensional form is