Diffusiophoresis of a Weakly Charged Dielectric Fluid Droplet in a Cylindrical Pore

1. Introduction

When a solute concentration gradient is established somehow in a solution system, the solutes will migrate downward this concentration gradient as predicted by the famous Fick’s law [1], whether the solutes are electrolytes or non-electrolytes. A corresponding osmosis flow is thus generated in the solution referred to as the diffusioosmosis flow [2]. When this flow runs into a colloidal entity suspended in the solution, it will set the colloidal entity in motion, referred to as the diffusiophoretic motion or diffusiophoresis, a very important phoretic motion driven by solute concentration gradient nearby. As the colloidal entities are generally charged and suspended in electrolyte solution, it is the ionic solutes in an electrolyte solution that are of major interest in practice, hence is chosen as the scope of the current study.

Compared with the well-known electrophoresis where the colloidal entity is driven by an external electric field applied upon the electrolyte solution, the diffusiophoresis is a lesser-known electrokinetic phenomenon but has been gaining increasing attention in recent years due to its unique features such as there is no or nearly negligible Joule heating effect in diffusiophoresis, for instance. This makes it extremely suitable in biomedical applications like drug delivery in human body, as an increase of four degrees Celsius is fatal to mammal cells in general. And it would be impractical, if not impossible to apply an external electric field upon the human body to manipulate the motion of the nanomedicines toward the desired region as it is often not known in advance and the human body is too complicated hence there is no way to follow the precise route of the nanomedicine there. Diffusiophoresis, on the other hand, is capable of reaching the region needing therapy on its own. This is because the injured or sick tissue in the human body often releases specific chemicals into the body fluid and establishes a local solute concentration gradient nearby, which attracts and guides the nanomachines to its neighborhood like a cruise missile following the guiding signals to its desired area. Moreover, diffusiophoresis has been utilized in enhanced oil recovery (EOR) as well, where the remaining crude oil droplets in the porous reservoir are driven out successfully via the solute concentration gradient induced by the injection of seawater, the electrolyte solution of NaCl. Moreover, many novel applications have been reported in recent years utilizing the diffusiophoresis mechanism, such as using diffusiophoresis as an electrokinetic mechanism for particle separations in microfluidic and nanofluidic operations [3,4,5,6,7].

The nanomedicines in drug delivery are often in the form of liposomes, where the therapeutic medicines are dissolved in the fluid filling the droplet with a lipid bilayer as the droplet surface. The migration of crude oil droplets in EOR certainly is a perfect example of droplet diffusiophoresis [3]. Note that a rigid particle and a gas bubble are just limiting examples of droplets with extremely high or low viscosities, respectively. Thus, the exploration of droplet phoretic motions will provide insights and information valid for rigid particles and gas bubbles as well. On the other hand, the cylindrical pore is widely used in conventional capillary electrophoresis in DNA sequencing [8,9,10,11], protein analysis[12,13,14,15], and so on, due to its convenience in operations and efficiency in temperature control to reduce the impact of Joule heating effect, among other things. As a matter of fact, phoretic motions of droplets have been frequently encountered in microfluidic and nanofluidic operations, due to their various merits [16,17,18,19,20]. The internal flow within the droplets helps mix the drugs efficiently, which can effectively reduce the cost of drug formulation. As a result, we decide to conduct a theoretical study on the diffusiophoretic motion of a dielectric droplet in a cylindrical pore in particular to understand the underlying electrokinetic mechanisms leading to the ultimate droplet motion, and find out the key parameters one can use to manipulate the droplet motion as desired. In other words, we extend the exploration to investigate the corresponding capillary diffusiophoresis. This classic geometric configuration certainly indicates its direct applications in micro/nanofluidic operations. As for the droplet of particular interest in drug delivery, the liposomes are often in the size range from 20 nanometers to 50 nanometers in practice, or even larger sometimes [21,22]. On the other hand, capillary with a diameter of 0.2 micron are commercially available nowadays, and pore size as small as 5 nanometers can be fabricated in lab [23,24], Hence the boundary confinement effect, both hydrodynamic and electrokinetic, due to the presence of a nearby channel wall in a very narrow pore has to be reckoned with carefully. As a result, this will be the focus of this study to investigate its impact upon the droplet motion, among other parameters of electrokinetic interest, such as the droplet size, the electrolyte strength, the radius of the cylindrical pore, and so on. In addition, it should be noted that a colloid in the cylindrical pore has been used as a modelling way for the migration of a colloid through a porous medium with the virtual pore size evaluated from the porosity of the medium [25]. Thus, the findings in the current study are applicable to the diffusiophoresis in a porous medium as well.

Phoretic motions of a colloidal entity in a cylindrical pore have been investigated theoretically by many research groups [26,27,28,29,30,31,32,33,34,35,36]. Most of them focus on the electrophoresis phenomenon though. The presence of a nearby cylindrical pore has been found to impose a profound boundary confinement effect upon the colloidal motion in general. The hydrodynamic drag force is increased due to the no-slip boundary condition of the cylindrical pore. Large viscous shear force is thus generated, which tends to slow down the droplet motion, for one thing. Moreover, the distribution of ions within the electric double layer, the source of the electric driving force setting the colloid in motion, can be severely altered as well if the cylindrical pore is sufficiently narrow. This in turn changes the electrostatic environment surrounding the droplet, hence the corresponding electric driving force. The narrower the cylindrical pore is, the more significant the boundary confinement is in general, both electrostatically and hydrodynamically. In particular, Lee and his coworkers [26] explored the electrophoresis of a weakly charged dielectric droplet in a cylindrical pore recently, which is the first time the phoretic motion of a dielectric droplet has ever been studied in a cylindrical pore. Interesting phenomena like mobility reversal were reported there. For droplet diffusiophoresis in a cylindrical pore, however, there has been no report in the literature, to the best of our knowledge.

Note in particular that the shear stress upon the droplet surface generated when a droplet is in motion is very detrimental to the integrity of the droplet. Because the droplet surface may be destroyed if the shear stress is beyond some critical magnitude that the surface can sustain. It has been reported that in highly constrict blood vessels, shear stress can increase by one to two orders of magnitude locally compared to normal vessels [37,38]. This may lead to the premature collapse of a liposome in drug delivery and loss of the encapsulated therapeutic medicines before it reaches the desired region needing therapy [39,40,41]. In addition, it may cause fatal breakups of the cell membranes for living cells, such as stem cells, which conduct phoretic motions in designed separation operations or in biomedical analyses in a lab-on-chip device [42,43]. As a result, it is highly desirable to operate the necessary maneuvers far away from the situation where large shear stress upon the surface of droplets may occur, as mentioned above. It turns out that the solidification phenomenon provides a promising design guideline in the fabrication stage of droplet nanomedicines like liposomes, as it will generate ideal and resilient drug carriers to achieve the maximum therapeutic performance in the human body. The details will be elaborated in the subsequent Results and Discussion section.

In summary, we theoretically investigate the diffusiophoretic motion of a dielectric droplet in a cylindrical pore in this study. Parameters of electrokinetic interest, like the droplet size, the electrolyte strength of the suspending solution, the ratio of the droplet viscosity to that of the ambient solution, and the ratio of the droplet radius to that of the cylindrical pore, are explored to investigate their respective effect and joint impact upon the droplet motion. In order to highlight the strong boundary confinement effect in a very narrow pore, we focus on the diffusiophoretic motion of a weakly charged dielectric droplet within a chargeless cylindrical pore here.

2. Theory

As shown in Figure 1, we consider the diffusiophoretic motion of a dielectric fluid droplet in response to a solute concentration gradient, ∇C, along the axis of a cylindrical pore. The cylindrical pore is filled with a binary electrolyte solution, such as KCl. The droplet is filled with chargeless dielectric fluid with a uniform surface charge distribution, assuming constant surface charge density, which remains constant with all the varying electrokinetic parameters, such as κa, where κ is the electrokinetic strength and a is the droplet radius. This is one of the classic surface charge conditions in addition to the constant surface potential and the charge-regulation condition [2,44]. Spherical coordinates (r, θ, φ) are adopted with the origin located at the center of the moving droplet for the domain inside the droplet. Cylindrical coordinates (R, Θ, Z), on the other hand, are adopted for the domain between the droplet surface and the cylindrical pore. The droplet moves along the axis of the cylindrical pore accordingly with a constant velocity U.

The electrolyte solution is assumed to be Newtonian and incompressible with constant viscosity and electric permittivity. The droplet surface is assumed ion-impenetrable. And it maintains its spherical shape as it moves, which is justified by the extremely low hydrodynamic Weber number, typically around 10⁻⁷ [45].

2.1. Governing Electrokinetic Equations

The governing equations are based on the general fundamental electrokinetic theory [44] which consists of the Poisson equation, the ions migration equations for both cations and anions (Nernst-Planck equation), the momentum equations for fluids both inside and outside the droplet, and the ions conservation equation for both cations and anions as shown below:

$${\nabla }^2\mathsf{\varphi }=-{\displaystyle \frac{\mathsf{\rho }}{{\mathsf{\epsilon }}_{\mathrm{m}}}}$$

(1)

$${\nabla }^2{\mathsf{\varphi }}_{\mathrm{I}}=-{\displaystyle \frac{\mathsf{\rho }}{{\mathsf{\epsilon }}_{\mathrm{D}}}}$$

(2)

$${\mathbf{f}}_{\mathbf{j}}=-{\mathrm{D}}_{\mathrm{j}}\left(\nabla {\mathrm{n}}_{\mathrm{j}}+{\displaystyle \frac{{\mathrm{z}}_{\mathrm{j}}\mathrm{e}}{{\mathrm{k}}_{\mathrm{B}}\mathrm{T}}}{\mathrm{n}}_{\mathrm{j}}\nabla \mathsf{\varphi }\right)+{\mathrm{n}}_{\mathrm{j}}\mathbf{v}$$

(3)

$$\nabla \cdot {\mathbf{f}}_{\mathbf{j}}=0$$

(4)

$${\mathsf{\eta }}_{\mathrm{m}}{\nabla }^2\mathbf{v}-\nabla \mathrm{P}-\mathsf{\rho }\nabla \mathsf{\varphi }=0$$

(5)

$${{\mathsf{\eta }}_{\mathrm{D}}\nabla }^2{\mathbf{v}}_{\mathrm{I}}-\nabla {\mathrm{P}}_{\mathrm{I}}=0$$

(6)

$$\nabla \cdot \mathbf{v}=0$$

(7)

$$\nabla \cdot {\mathbf{v}}_{\mathrm{I}}=0$$

(8)

where eqns. (1) and (2) are the Poisson equations based on Gauss’s divergence theorem for fluids both interior and exterior to the droplet, with ϕ denoting the local electrical potential and $\mathsf{\rho }=\sum _{\mathrm{j}=1}^{\mathrm{N}}{\mathrm{z}}_{\mathrm{j}}\mathrm{e}{\mathrm{n}}_{\mathrm{j}}$ the space charge density, the total number of electric charges per unit volume. N is the number of the ion species in the fluid. Moreover, n_j refers to the number concentration of ion species j, z_j is its valence number, and D_j its diffusivity coefficient. Eq. (3) is the ion migration equation for each ion species j, the Nernst–Planck equation. Equations (5) and (6) are the momentum equations of the fluid’s exterior and interior to the droplet, respectively, with P standing for the pressure, η_m and η_D the viscosities of the ambient solution and interior fluid, respectively. Eqns. (7) and (8) are the incompressibility constraints for fluids both inside and outside the droplet. In Eq. (5), an additional term representing the electric body force in the surrounding electrolyte solution, -ρ∇ϕ, is included, which is the physical origin of the electric Maxwell stress tensor responsible for the local electric force. R_w is the radius of the cylindrical pore. In addition, εm is the electric permittivity of the ambient electrolyte solution, whereas ε_D is the electric permittivity of the droplet fluid. Suffix I overall indicates the region inside the droplet. The definitions of the rest of the symbols can be found in the List of Symbols in the Supplementary Material. Note that the above governing equations are written in vector form, hence are independent of the specific coordinates adopted.

Before the imposition of the solute concentration gradient, the system is at equilibrium, and a set of equilibrium solutions is obtained based on the above treatment. Once the concentration gradient is applied, the system is disturbed and all the dependent variables will deviate from their equilibrium values. Assuming the concentration gradient is very small compared to the equilibrium concentration, these deviations will be small as well. Standard linear perturbation analysis thus can be adopted [2,46]. Assuming the variables after the disturbance satisfy the above governing equations as well, and subtracting the equilibrium equations from them, an extra set of equations governing the perturbations of the variables is obtained.

The complete governing equations are as follows:

$${\nabla }^{\mathrm{*}2}{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}+{\displaystyle \frac{{\left(\mathsf{\kappa }\mathrm{a}\right)}^2}{1+\mathsf{\alpha }}}[\exp \left(-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}\right)-\exp \left({\mathsf{\alpha }\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}\right)]=0$$

2.2. Boundary Conditions

The following are the associated boundary conditions on the droplet surface and at infinity, where the concentration gradient is imposed upon the system:

$${\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}=-{\mathsf{\sigma }}^{\mathrm{*}}$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(17)

$${\left.{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\right|}_{\mathrm{r}\mathrm{*}=1}=0$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(18)

$${\displaystyle \frac{\partial {\mathrm{g}}_1^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}=0$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(19)

$${\displaystyle \frac{\partial {\mathrm{g}}_2^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}=0$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(20)

$${\left.{\displaystyle \frac{\partial {\mathsf{\Psi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\right|}_{\mathrm{r}\mathrm{*}=1^{+}}={\left.{\displaystyle \frac{\partial {\mathsf{\Psi }}_{\mathrm{I}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\right|}_{\mathrm{r}\mathrm{*}=1^{-}}$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(21)

where Eqn. (17) indicates that the dielectric droplet under consideration here possesses constant surface charge density. Eqn. (18), on the other hand, indicates that the droplet is dielectric which reduces to the electric insulating condition, i.e., the right-hand side is set to zero, if the electric permittivity of the interior fluid is much smaller than that of the ambient electrolyte solution, such as the silicone oil droplet in an aqueous electrolyte solution assumed here for simplicity. Eqns. (19) and (20) are derived based on the assumption that the droplet is ion-impenetrable. Eqn. (21) indicates that the tangent fluid velocities across the droplet surface should be continuous.

$${{(\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }}^{\mathbf{N}}+{{\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }}^{\mathbf{M}}{\left)\right|}_{{\mathrm{r}}^{\mathrm{*}}=1^{+}}={({\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }\mathrm{I}}^{\mathbf{N}}+{{\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }\mathrm{I}}^{\mathbf{M}}){|}_{{\mathrm{r}}^{\mathrm{*}}=1^{-}}$$

$${\mathrm{r}}^{\mathrm{*}}=1$$

(22)

Eqn. (22) is the generalized Rybczynski-Hadamard condition taking into account the electric Maxwell shear stress on the droplet surface from the exterior electrolyte solution [47]. It states that the total shear stress, including both hydrodynamic shear stress,${{\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }}^{\mathbf{N}}$, and electrostatic Maxwell shear stress, ${{\mathsf{\tau }}^{\mathrm{*}}}_{\mathrm{r}\mathsf{\theta }}^{\mathbf{M}}$, should be continuous across the droplet surface in the absence of the interfacial tension [47,48,49]. Complete Maxwell stress tensor in general is defined as follows: ${\mathsf{\tau }}^{\mathbf{M}}=\mathsf{\epsilon }\left[\nabla \mathsf{\varphi }\nabla \mathsf{\varphi }-{\displaystyle \frac{1}{2}}\left(\nabla \mathsf{\varphi }\bullet \nabla \mathsf{\varphi }\right)\mathbf{I}\right]$, where I is the identity tensor.

The boundary conditions at the cylindrical pore surface are listed below in dimensionless form:

$${\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}={\mathsf{\varphi }}_{\mathrm{w}}^{\mathrm{*}}=0$$

(23)

$${\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{R}}^{\mathrm{*}}}}=0$$

(24)

$${\displaystyle \frac{\partial {\mathrm{g}}_{\mathbf{j}}^{\mathrm{*}}}{\partial {\mathrm{R}}^{\mathrm{*}}}}=0$$

(25)

$${\mathsf{\psi }}^{\mathrm{*}}=0$$

(26)

$${\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial {\mathrm{R}}^{\mathrm{*}}}}=0$$

(27)

Moreover, far away from the droplet, both upstream and downstream, the flow field and the electric field should approach asymptotically the corresponding situation in the absence of the droplet as shown below:

$${\mathsf{\psi }}^{\mathrm{*}}={\displaystyle \frac{1}{2}}{\mathrm{U}}^{\mathrm{*}}{\mathrm{R}}^{\mathrm{*}2}+{\mathsf{\psi }}_{\mathrm{\infty }}^{\mathrm{*}}\left({\mathrm{R}}^{\mathrm{*}}\right),{\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial {\mathrm{Z}}^{\mathrm{*}}}}=0$$

$${\mathbf{Z}}^{\mathrm{*}}=\pm {\mathbf{L}}^{\mathrm{*}}$$

(28)

$${\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}={\mathsf{\varphi }}_{\mathrm{e},\mathrm{\infty }}^{\mathrm{*}}\left({\mathrm{R}}^{\mathrm{*}}\right)=0$$

$${\mathbf{Z}}^{\mathrm{*}}=\pm {\mathbf{L}}^{\mathrm{*}}$$

(29)

$${\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{Z}}^{\mathrm{*}}}}=-\mathsf{\beta }{\nabla }^{\mathrm{*}}{\mathrm{C}}^{\mathrm{*}}=0$$

$${\mathbf{Z}}^{\mathrm{*}}=\pm {\mathbf{L}}^{\mathrm{*}}$$

(30)

$${\mathrm{g}}_{\mathbf{j}}^{\mathrm{*}}=\left(\mathsf{\beta }-1\right){\nabla }^{\mathrm{*}}{\mathrm{C}}^{\mathrm{*}}{\mathrm{R}}^{\mathrm{*}}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }=-{\nabla }^{\mathrm{*}}{\mathrm{C}}^{\mathrm{*}}{\mathrm{R}}^{\mathrm{*}}\mathrm{c}\mathrm{o}\mathrm{s}\mathsf{\theta }$$

$${\mathbf{Z}}^{\mathrm{*}}=\pm {\mathbf{L}}^{\mathrm{*}}$$

(31)

where L^*(=L/a) stands for the scaled distance large enough so that the droplet mobility does not change with increasing L^* anymore. It turns out that L^* = 10 is sufficient. A dimensionless index β, defined as β $={\displaystyle \frac{{\mathrm{D}}_1-{\mathrm{D}}_2}{{\mathrm{D}}_1+{\mathrm{D}}_2}}$ in a binary electrolyte solution appears in Eqn. (30), where D₁ is the diffusivity of the cations and D₂ is the diffusivity of the anions. β is a measurement of the strength of the induced diffusion potential [50,51]. For a KCl solution with approximately identical diffusivities of cations and anions, β is often regarded as zero as a benchmark situation of diffusion potential, which is set to zero here as well. Eqn. (31) indicates the impact of the diffusion potential upon the droplet motion, which is derived based on the electro-neutrality constraint in diffusiophoresis. Here, β is set to zero in Eqn. (30) and (31) in a KCl solution. Eqns. (30) and (31) are the boundary conditions at infinity in general. They indicate that there should be no net electric current at infinity since there is no external electric field applied upon the system

2.3. Mathematical Treatments

As the system is axisymmetric, there is no dependence on φ and the problem reduces to a two-dimensional one. Moreover, as the droplet interior fluid is chargeless, there is no need to solve for the electric field there. Since there is an analytical formula for the purely hydrodynamic fluid field inside a droplet, which can be incorporated into the corresponding exterior fluid field as a boundary condition, only the electric and fluid flow fields in the exterior region need to be solved. To accomplish this, a conformal mapping technique is adopted to convert the irregular domain exterior to the droplet into three consecutive rectangular sub-domains with continuity conditions of all the system variables and their derivatives assumed across the boundaries between each of these sub-domains. A patched pseudo-spectral method based on Chebyshev polynomials is then applied in both of the orthogonal directions of each subdomain [52,53,54]. Detailed mathematical treatments can be found elsewhere [2,55]. A numerical algorithm is then developed based on the above treatments, which is used to solve for the coupled electric and flow fields of interest here.

2.4. Evaluation of Droplet Mobility

Once the coupled electric field and the flow field are solved via the numerical algorithm developed, the corresponding dimensionless hydrodynamic drag force (${\mathrm{F}}_{\mathbf{D}\mathbf{z}}^{\mathrm{*}}$) and the electric driving force ${\mathrm{F}}_{\mathbf{E}\mathbf{z}}^{\mathrm{*}}$ upon the droplet surface can be evaluated, respectively, as follows:

$${\mathrm{F}}_{\mathbf{D}\mathbf{z}}^{\mathrm{*}}=\mathsf{\pi }{\int }_0^{\mathsf{\pi }}{\left|{\mathrm{r}}^{\mathrm{*}4}{\sin }^3\mathsf{\theta }{\displaystyle \frac{\partial }{\partial {\mathrm{r}}^{\mathrm{*}}}}\left({\displaystyle \frac{{\mathrm{E}}^{\mathrm{*}2}{\mathsf{\psi }}^{\mathrm{*}}}{{\mathrm{r}}^{\mathrm{*}2}{\sin }^2\mathsf{\theta }}}\right)\right|}_{{\mathrm{r}}^{\mathrm{*}}=1}\mathrm{d}\mathsf{\theta }-\mathsf{\pi }{\int }_0^{\mathsf{\pi }}{\left|{\mathrm{r}}^{\mathrm{*}2}{\sin }^2\mathsf{\theta }\left\{{\displaystyle \frac{(\mathsf{\kappa }\mathrm{a}{)}^2}{1+\mathsf{\alpha }}}\left[\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})\left(1-\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}-{\mathrm{g}}_1^{\mathrm{*}}\right)-\exp (\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})(1+\mathsf{\alpha }(\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}+{\mathrm{g}}_2^{\mathrm{*}}\left)\right)\right]{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}+\left[\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})-\exp (\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})\right]{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial \mathsf{\theta }}}\right\}\right|}_{{\mathrm{r}}^{\mathrm{*}}=1}\mathrm{d}\mathsf{\theta }$$

(32)

$${\mathrm{F}}_{\mathbf{E}\mathbf{z}}^{\mathrm{*}}=\mathsf{\pi }{\int }_0^{\mathsf{\pi }}{\left|\left\{{\mathrm{r}}^{\mathrm{*}2}\sin \mathsf{\theta }\cos \mathsf{\theta }\left[{\left({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\right)}^2+2{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\right]+{\mathrm{r}}^{\mathrm{*}2}{\sin }^2\mathsf{\theta }\left({\displaystyle \frac{{\partial }^2{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}\partial \mathsf{\theta }}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}+{\displaystyle \frac{{\partial }^2{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}\partial \mathsf{\theta }}}{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}+{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{{\partial }^2\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}\partial \mathsf{\theta }}}+{\displaystyle \frac{1}{{\mathrm{r}}^{\mathrm{*}2}}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{{\partial }^2{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathsf{\theta }}^2}}+{\displaystyle \frac{1}{{\mathrm{r}}^{\mathrm{*}2}}}{\displaystyle \frac{{\partial }^2{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathsf{\theta }}^2}}{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial \mathsf{\theta }}}+{\displaystyle \frac{1}{{\mathrm{r}}^{\mathrm{*}2}}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{{\partial }^2\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathsf{\theta }}^2}}\right)-2{\mathrm{r}}^{\mathrm{*}}{\sin }^2\mathsf{\theta }\left({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}+{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial \mathsf{\theta }}}+{\displaystyle \frac{\partial \mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}\right)\right\}\right|}_{{\mathrm{r}}^{\mathrm{*}}=1}\mathrm{d}\mathsf{\theta }$$

(33)

The droplet diffusiophoretic mobility, defined as droplet velocity divided by the magnitude of the concentration gradient imposed, can be evaluated as follows:

$${\mathsf{\mu }}^{\mathrm{*}}={\displaystyle \frac{{\mathbf{U}}^{\mathrm{*}}}{{\nabla }^{\mathrm{*}}{\mathrm{C}}^{\mathrm{*}}}}=-{\displaystyle \frac{{\mathbf{F}}_{\mathbf{E}}^{\mathrm{*}}(\mathrm{o}\mathrm{b}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{u}\mathrm{b}-\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{t}\mathrm{w}\mathrm{o})}{{\mathbf{F}}_{\mathbf{D}}^{\mathrm{*}}(\mathrm{o}\mathrm{b}\mathrm{t}\mathrm{a}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{n}\mathrm{s}\mathrm{u}\mathrm{b}-\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{n}\mathrm{e})}}$$

(34)

where the ingenious methodology by O’Brien and White [46] is adopted which demonstrates that the system can be decomposed into two separate artificial auxiliary sub-problems, with one of them corresponding to a purely electrostatic one, sub-problem one, and the other purely hydrodynamic one, sub-problem two.

$${\nabla }^{\mathrm{*}2}{\mathsf{\varphi }}_{{\mathrm{e}}_{\mathrm{I}}}^{\mathrm{*}}=0$$

$${{\nabla }^{\mathrm{*}}}^2\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}-{\displaystyle \frac{(\mathsf{\kappa }\mathrm{a}{)}^2}{1+\mathsf{\alpha }}}\left[\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})+\mathsf{\alpha exp}(\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})\right]\mathsf{\delta }{\mathsf{\varphi }}^{\mathrm{*}}={\displaystyle \frac{(\mathsf{\kappa }\mathrm{a}{)}^2}{1+\mathsf{\alpha }}}\left[\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}){\mathrm{g}}_1^{\mathrm{*}}+\mathsf{\alpha }\exp (\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}){\mathrm{g}}_2^{\mathrm{*}}\right]$$

$${\nabla }^{\mathrm{*}2}\mathsf{\delta }{\mathsf{\varphi }}_{\mathrm{I}}^{\mathrm{*}}=0$$

where$\mathsf{\delta }$refers to a perturbation amount of the system variable after it. Eqn. (9) is the governing equation for the equilibrium electric potential in the ambient electrolyte solution. Eqn. (10), on the other hand, indicates that there are no electric charges or ions in the interior droplet fluid. ${\mathrm{R}}_{\mathrm{w}}^{\mathrm{*}}$is the ratio of the cylindrical pore radius to the droplet radius, i.e., R_w/a. Superscript * indicates dimensionless physical quantities here in general, with specific definitions contained in the List of Symbols in the Supplementary material. Eqn. (11) and Eqn. (12) are the corresponding governing equations of the electric potential disturbance based on the standard linear perturbation analysis.$\mathsf{\alpha }$is unity for a symmetric binary electrolyte solution considered here. g_j is the shape function showing deviation from the concentric Boltzmann distribution of ion species j when the droplet is in motion due to the effect of convection flow alone. Precise definitions can be found in the List of Symbols. Detailed further elaborations can be found elsewhere [2]. Hence

$${{\nabla }^{\mathrm{*}}}^2{\mathrm{g}}_1^{\mathrm{*}}-({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathrm{g}}_1^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}+{\displaystyle \frac{1}{{{\mathrm{r}}^{\mathrm{*}}}^2}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{\partial {\mathrm{g}}_1^{\mathrm{*}}}{\partial \mathsf{\theta }}})={\displaystyle \frac{{\mathrm{P}\mathrm{e}}_1}{{{\mathrm{r}}^{\mathrm{*}}}^2\sin \mathsf{\theta }}}({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}-{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial \mathsf{\theta }}})$$

$${{\nabla }^{\mathrm{*}}}^2{\mathrm{g}}_2^{\mathrm{*}}+\mathsf{\alpha }({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathrm{g}}_2^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}+{\displaystyle \frac{1}{{{\mathrm{r}}^{\mathrm{*}}}^2}}{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{\partial {\mathrm{g}}_2^{\mathrm{*}}}{\partial \mathsf{\theta }}})={\displaystyle \frac{{\mathrm{P}\mathrm{e}}_2}{{{\mathrm{r}}^{\mathrm{*}}}^2\sin \mathsf{\theta }}}({\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}{\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}-{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}{\displaystyle \frac{\partial {\mathsf{\psi }}^{\mathrm{*}}}{\partial \mathsf{\theta }}})$$

where Eqns. (13) and (14) are the corresponding Nernst-Planck equations governing the migration of cations and anions, respectively, in response to the solute concentration gradient. In addition,

$${\mathrm{E}}^4{\mathsf{\psi }}^{\mathrm{*}}=-{\displaystyle \frac{(\mathsf{\kappa }\mathrm{a}{)}^2}{1+\mathsf{\alpha }}}\left\{\begin{array}{c}\left[{\displaystyle \frac{\partial {\mathrm{g}}_1^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})+{\displaystyle \frac{\partial {\mathrm{g}}_2^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\mathsf{\alpha }\exp (\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})\right]{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial \mathsf{\theta }}}\\ -\left[{\displaystyle \frac{\partial {\mathrm{g}}_1^{\mathrm{*}}}{\partial \mathsf{\theta }}}\exp (-{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})+{\displaystyle \frac{\partial {\mathrm{g}}_2^{\mathrm{*}}}{\partial \mathsf{\theta }}}\mathsf{\alpha }\exp (\mathsf{\alpha }{\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}})\right]{\displaystyle \frac{\partial {\mathsf{\varphi }}_{\mathrm{e}}^{\mathrm{*}}}{\partial {\mathrm{r}}^{\mathrm{*}}}}\end{array}\right\}\sin \mathsf{\theta }$$

$${\mathrm{E}}^4{\mathsf{\psi }}_{\mathrm{I}}^{\mathrm{*}}=0$$

where Eqns. (15) and (16) are the corresponding momentum equations governing the flow field outside and inside the droplet, respectively. Stream functions in dimensionless form, ${\mathsf{\psi }}^{\mathbf{*}}$ and ${\mathsf{\psi }}_{\mathrm{I}}^{\mathrm{*}}$ are introduced to eliminate the pressure term and satisfy the incompressible constraints, Eqns. (7) and (8), automatically.

3. Results and Discussion

We first check the convergence of the numerical scheme in that the calculation results do not change any more with further mesh refinement. It turns out that 51 grid points in the θ-direction plus 100 grid points in the r-direction before conformal mapping are sufficient. We also check the convergence behavior of the droplet mobility with increasing radius of the cylindrical pore. Indeed, they all approach asymptotically the single droplet situation in the absence of the cylindrical pore, as shown in Figure 2, where the droplet mobilities as functions of the ratio of the droplet radius to the cylindrical pore radius are shown for several viscosities of the droplet interior fluids. We thus conclude that the numerical scheme is reliable and the calculation results are accurate. We then go on to investigate the droplet motion based on it. Note that, as shown in Figure 2, mobility reversal is observed for a highly viscous droplet, which is essentially a rigid particle with the ratio of the droplet viscosity to that of the ambient solution equal to 100. This indicates that the moving direction of a droplet in narrow channels may be opposite to a droplet in wide channels, which has potential applications for droplet manipulation in microfluidic and nanofluidic devices.

at various viscosity ratio σ_H for a dielectric droplet with σ*=2.03 and κa=1 in the KCl solution (β=0): (A) small scale; (B) large scale.

Figure 3 demonstrates the mobility profiles of a benchmark weakly charged droplet with the dimensionless surface charge density σ^* equals to 2.03 in a very narrow cylindrical pore with ${\mathrm{R}}_{\mathrm{w}}^{\mathrm{*}}$ =1.2 as functions of κa for various viscosity ratios of the droplet interior fluid to that of the ambient electrolyte solution, where κ is the electrolyte strength and a is the droplet radius. κa can be regarded as a measurement of the electric double layer thickness. The larger the value of κa is, the thinner the electric double layer surrounding the droplet and vice versa. This is due to the double layer compression effect when the electrolyte strength κ is high in the bulk solution. In fact, κa can be regarded as the dimensionless reciprocal of the double layer thickness, since κ^-1, the Debye length, is a characteristic length of the double layer [44]. It is interesting to note that for a fixed droplet size in such a narrow channel, the droplet tends to move upward against the concentration gradient in relatively dilute electrolyte solutions, whereas downward along the concentration gradient in relatively concentrated solutions. The underlying fundamental electrokinetic mechanisms are elaborated as follows.

= 1.2 with σ^*=2.03 in the KCl solution (β=0).

When a concentration gradient of ions is somehow established in an electrolyte solution, the ions will migrate downward this concentration gradient in the bulk solution by the diffusion mechanism instantly, both the cations and the anions. An osmosis flow referred to as the diffusioosmosis flow is thus generated due to the closely bonded hydrated solvent molecules moving together with the ions and the momentum transfer to the pure solvent molecules in the neighborhood via the hydrodynamic stress, as indicated by the viscosity of the solution. If the diffusivities of the cations and anions are the same, they will migrate at the same speeds. If they are different, an induced diffusion potential will be instantly generated within the bulk electrolyte solution by the Coulomb electrostatic force, which tends to speed up the slower ions and slow down the faster ones. Hence eventually, and almost instantly, they move downward the concentration gradient in the bulk electrolyte solution at the same speeds as well. The diffusiophoresis observed in the former electrolyte solution is called chemiphoresis sometimes, emphasizing the fact that there is no induced diffusion potential involved and the phoretic motion of the colloidal entity is purely due to the concentration/chemical affinity of the ions. In the latter case, in addition to the motion caused by the chemiphoresis mechanism, the induced diffusion potential inside the bulk electrolyte solution tends to behave as an inner battery to drive the charged colloidal entity just like an externally applied electric field. The motion induced by this mechanism is called the electrophoresis component of the diffusiophoresis, emphasizing the electrokinetic origin of its driving mechanism. Note that the chemophoresis mechanism is always present with or without the electrophoresis mechanism. It is coupled with the electrophoresis component if there is indeed a diffusion potential generated in the electrolyte solution. Hence it is sometimes referred to as the chemiphoresis component. Note that these two components are coupled and cannot be rigorously separated from each other. When these ions in the diffusioosmosis flow run into a charged droplet ahead, both cations and the anions migrate across the boundary of the electric double layer surrounding the charged droplet at the same speed and total amount. The counterions will be predominantly attracted in the upstream hemisphere of the charged droplet due to the Coulomb electrostatic attraction force between them and the droplet surface charges. The coions, on the other hand, will be predominantly repulsed and swept toward the downstream hemisphere instead, based on the Coulomb electrostatic repulsion force. An induced asymmetric double layer polarization is results which generates a local electric field driving the droplet in motion. The force balance between this electric driving force and the hydrodynamic drag force determines the ultimate droplet motion pattern. Note that while the fundamental scenario presented above indicates initially this electric driving force tends to drive the droplet upward against the concentration gradient based on the Coulomb electrostatic force, the actual moving direction of the droplet can go either way. Rigorous solution of the governing fundamental electrokinetic equations is the only way to sort it out, as presented in Figure 3 above.

When the ions in the downward-moving diffusioosmosis flow are relatively dilute, its associated downward pushing force is relatively weak as well. The upward electric driving force is hence dominant and pushes the droplet upward. A positive droplet mobility is thus observed, as shown in Figure 3. When κa increases further, however, the downward osmosis convection flow becomes significantly enhanced. The droplet is pushed downward eventually, and a negative mobility is observed as shown in Figure 3 as well. This indicates an interesting phenomenon that one can manipulate the moving direction of a droplet in a narrow channel simply by changing the concentration of the electrolyte solution. Moreover, under the same geometric configuration filled with identical electrolyte solution, smaller droplets tend to move to the region with higher concentration of ions and larger ones in the opposite direction. In drug delivery with liposome droplets, the injured tissues or regions needing therapy often release specific chemicals in their neighborhood. A concentration gradient of these chemicals is thus established, which serves as the natural guiding and driving mechanisms to lead the liposomes to their desired final destination [56,57]. Smaller liposomes would be favorable according to Figure 3.

The boundary confinement effect is clearly demonstrated by direct comparison with the corresponding mobility profiles of a droplet suspended in an infinite medium of electrolyte solution, as shown in Figure 4. The droplet mobility is significantly increased at smaller κa, or thicker double layer, due to the presence of a nearby cylindrical pore. The ion concentration in the double layer is significantly increased by this geometric confinement effect, leading to a much stronger motion-inducing double layer polarization. The electric driving force is thus significantly enhanced, leading to a mobility behavior completely opposite to the single droplet situation. When the double layer is thin, indicated by a relatively large κa value, however, the nearby channel wall will not be able to confine the double layer as profoundly as before. The boundary confinement effect will weaken significantly. The droplet mobility profiles, hence, assume a pattern similar to the corresponding single droplet situation.

= 1.2 in the KCl solution (β=0): (A) diffusiophoresis; (B) electrophoresis.

Other than the obvious confinement effect upon the droplet mobility, the presence of the nearby channel wall also has a subtle impact on the ions distribution upstream and downstream of the droplet. Due to the significant reduction of area available for the downward diffusioosmosis flow at the equator area, the ions have a huge buildup in the upstream. This in turn induces an extra outer region of a secondary double layer based on the electrostatic Coulomb law. An extra enhancement of coions results downstream as well, which also induces an extra outer region enriched in counterions. Moreover, as the macro/convection osmosis flow has to pass through the narrow annulus region at the equator, severe buildup of ions also takes place near the channel wall, both upstream and downstream of the droplet. The above deduction based on electrokinetics is clearly observed in Figure 5(A), where the precise shape of the double layer polarization is demonstrated. Compared with the corresponding double layer polarization in droplet electrophoresis, as shown in Figure 5(B), the ion distribution is much more complicated in diffusiophoresis.

It is also interesting to note in Figure 3 and Figure 6(C), there is a critical value of κa where all droplets move at identical speeds regardless of their viscosities. They move like rigid particles in that the characteristic recirculating interior vortex flows pertinent to typical droplet motion completely disappear. We thus call it the “solidification phenomenon” [57,58]. This intriguing phenomenon is due to the deadlock of the spinning electric Maxwell traction and the hydrodynamic drag force on the droplet surface, which are in opposite directions here [57,58]. Further detailed electrokinetic mechanisms will be elaborated as follows by examining more droplet motion patterns, both droplets moving directions and surface spinning orientations.

= 1.2, in the KCl solution (β=0): (A) κa=0.5; (B) κa=1.55; (C) κa=1.69; (D) κa=2.

It is interesting to note that while the droplet moves upward, indicating the overall dominance of the upward electric driving force relative to the weak downward osmosis flow, the portion attributes to the spinning force on the droplet surface may be still no match to the hydrodynamic drag force spinning droplet surface downward, as shown in Figure 6(A). The major portion of the electric driving force goes to the pressure drag (form drag) pushing the droplet upward as a whole by exerting the force in the normal direction of the droplet surface [47]. As a result, the remaining viscous drag (skin drag), which spins the droplet surface in the tangent direction, is insufficient in this case to reverse the original spinning orientation on the droplet surface by hydrodynamic drag force from the osmosis flow [47]. However, as shown in Figure 6(B), when κa further increases to 1.4 the droplet moves downward with an axisymmetric exterior vortex flow to reconcile the two electric and hydrodynamic shear forces spinning the droplet surface in opposite directions [58]. It indicates that the orientation of the surface spinning direction cannot be decided by the moving direction of the droplet alone. They are not correlated in a simple way in general. It can only be determined by the precise relative portions of the electric driving force exerted upon the droplet surface in the form of either the pressure drag or the viscous drag upon the droplet. This distribution of relative portions is case-dependent and can only be determined by solving the governing electrokinetic equations coupled with associated boundary conditions rigorously as presented here. There is no simple rule of thumb for the prediction of surface spinning direction and moving direction of a droplet as well as its speed for a complicated system, like the one considered here. As κa increases further, the opposite spinning forces from the electric driving force and hydrodynamic drag force become identical in magnitude eventually at a critical value of κa, and the droplet surface spins no more. Without the net surface spinning motion, the droplet now moves like a rigid particle without the interior recirculating vortex flow, which is referred to as the “solicitation phenomenon” here. As a result, the droplet mobility is independent of the droplet viscosity, as shown in Figure 6(C). In other words, regardless of the droplet’s viscosity, its mobility remains the same. Further increase of κa may even lead to the complete suppression of counterclockwise electric spinning motion via the overwhelming clockwise hydrodynamic spinning force, as shown in Figure 6(D).

The solidification phenomenon has important potential in practical applications like cell manipulation and drug delivery, to name a few. For animal cells and certain algal cells that lack a cell wall, for example, the cell contents are separated from the external environment only by the very thin lipid cell membrane. Cellular damage could be significantly enhanced by both the increase in intensity and the length of time over which the shear stress is applied [43]. Since the shear rate at the droplet surface is zero at the solidification point, the droplet can resist, or at least minimize, the undesirable external shear stress, thus avoiding potential cell damage during its motion. Another example can be found in typical bioreactor operations involving cells, where the vulnerability of cells to the shear stress generated by agitated mixing, among other sources, is determined by their resilience to it [42,43]. Additionally, in drug delivery with liposomes, the drug might leak during transportation because of the shear stress experienced before they reach the intended region needing therapy. This unwanted early release of the encapsulated nanomedicines can result in reduced overall drug performance and negative side effects [39,41,59]. With the solidification phenomenon pertinent to droplets, however, it provides the key information to solve all the problems mentioned above caused by the undesirable shear stress damaging to the droplet surface: the critical κa value at which the solidification phenomenon happens. For a cell of specific size, it provides a promising electrolyte strength at which the shear stress is absent, if diffusiophoresis is the driving mechanism for its motion. On the other hand, for drug delivery in human body, where the electrolyte strength is not an adjustable parameter, the optimal size of the drug-carrying liposomes in the fabrication stage can be obtained via the corresponding critical κa value, as diffusiophoresis is a crucial mechanism to enhance the concentration of drugs in the intended region of human body [60]. The reduction of cell mortality and enhancement of the success in drug delivery is anticipated in this way. Note that the solidification phenomenon is observed in the diffusiophoresis of a highly charged dielectric fluid droplet in an infinite medium of electrolyte solutions as well [45,57]. Here we find that it can happen for a weakly charged droplet in a narrow cylindrical pore as well. The enhancement of the motion-inducing double layer polarization by the boundary confinement effect in narrow channels is the underlying electrokinetic mechanism.

4. Conclusion

In this study, the diffusiophoretic motion of a weakly charged dielectric droplet in a cylindrical pore filled with KCl electrolyte solution is investigated theoretically. The droplet contains no electrolytes inside and has a constant, uniform surface charge density. Corresponding governing fundamental electrokinetic equations are solved with a patched pseudo-spectral method based on Chebyshev polynomials, coupled with a geometric mapping scheme to take care of the irregular shape of the solution domain. The results are summarized as follows:

Mobility reversal is observed in a benchmark narrow cylindrical pore. Moreover, an interesting and peculiar phenomenon referred to as the “solidification phenomenon” is observed at some specific critical droplet size or electrolyte strength in narrow channels under which the droplet mobilities are identical regardless of the droplet viscosities. The droplets move like rigid particles without the surface spinning motions or interior recirculation flows. As the corresponding shear rate is zero at this point, and so is the total shear stress upon the droplet surface, the droplet is resilient to exterior damaging shear stress tending to destroy the integrity of the droplet in motion. This provides an ideal and optimal droplet size in practical applications such as drug delivery, where the damage to the therapeutics-carrying liposomes can be minimized as a result before the nanomedicines reach the intended region of the human body. The overall performance can be maximized in this sense. For microfluidic and nanofluidic operations involving droplets, this provides a guideline for the ideal electrolyte strength of the solution filling the cylindrical pore as well, if the damage to the moving droplet by surface shear stress is of major concern, such as in cell operations. This solidification phenomenon is found in a highly charged dielectric droplet conducting either electrophoretic or diffusiophoretic motion in an infinite medium of electrolyte solution as well. We show here that with the presence of a very close channel wall, the same phenomenon can take place for a weakly charged droplet as well.

While the droplet may move either up against the concentration gradient or downward with it, there is no simple correlation between its moving direction and the orientation of its surface spinning direction. It can go either way. The relative distribution between the pressure drag and the viscous drag of the electric driving force is crucial to it. The force balance between the electric driving force, which is induced by the double layer polarization, and the hydrodynamic drag force from either side of the droplet surface determines the ultimate droplet motion pattern.

The precise shape of the motion-inducing double layer polarization is presented, which is much more complicated than the corresponding electrophoresis system. The convection downward diffusioosmosis flow is found to be responsible for it. Its impact is significantly enhanced in narrow channels in particular due to the speedup of the fluid flow through the narrow annular cross-section area surrounding the droplet equator. Buildup of ions upstream due to this reduction of fluid flow area takes place as a result. Induced secondary double layer polarization is observed under certain electrokinetic circumstances.

Size-dependent mobility in a cylindrical pore is observed, which provides a potential separation scheme in nano/microfluidic operations involving dielectric droplets of varying sizes.

Overall, the findings in this study provide a comprehensive understanding of the droplet diffusiophoresis behavior in a cylindrical pore, which is a classic geometric scheme in various practical applications involving droplets, such as drug delivery and droplet microfluidic/nanofluidic operations.