*Article* **Effect of Shear Flow on Nanoparticles Migration near Liquid Interfaces**

**Ali Daher 1,2,\*, Amine Ammar 1, Abbas Hijazi <sup>2</sup> and Lazhar Benyahia <sup>3</sup>**


**Abstract:** The effect of shear flow on spherical nanoparticles (NPs) migration near a liquid–liquid interface is studied by numerical simulation. We have implemented a compact model through which we use the diffuse interface method for modeling the two fluids and the molecular dynamics method for the simulation of the motion of NPs. Two different cases regarding the state of the two fluids when introducing the NPs are investigated. First, we introduce the NPs randomly into the medium of the two immiscible liquids that are already separated, and the interface is formed between them. For this case, it is shown that before applying any shear flow, 30% of NPs are driven to the interface under the effect of the drag force resulting from the composition gradient between the two fluids at the interface. However, this percentage is increased to reach 66% under the effect of shear defined by a Péclet number *Pe* = 0.316. In this study, different shear rates are investigated in addition to different shearing times, and we show that both factors have a crucial effect regarding the migration of the NPs toward the interfacial region. In particular, a small shear rate applied for a long time will have approximately the same effect as a greater shear rate applied for a shorter time. In the second studied case, we introduce the NPs into the mixture of two fluids that are already mixed and before phase separation so that the NPs are introduced into the homogenous medium of the two fluids. For this case, we show that in the absence of shear, almost all NPs migrate to the interface during phase separation, whereas shearing has a negative result, mainly because it affects the phase separation.

**Keywords:** liquid–liquid interface; shear rate; nanoparticles; diffuse interface; phase field method; molecular dynamics; numerical simulation

#### **1. Introduction**

In recent years, using nanoparticles (NPs) in industrial and medical markets has grown significantly. They have been of immense significance in different branches of science and engineering. This is basically due to their unique properties, such as augmented reactivity and special optical properties, which make them very suitable for products and applications in tissue engineering, composite technology, enhanced oil recovery and drug delivery [1,2]. In addition, NPs arise in nanoparticle-armored fluid droplets [3] and phase-arrested gels [4]. It is well known that studying biological processes on the nanoscale level is an essential point behind the development of nanotechnology [5].

The assembly of NPs at the liquid–liquid interface is essential in the preparation and stabilization of conventional emulsions, which have wide applications in petroleum, cosmetics, food, and biological transferring [6,7]. Modeling the dynamics of NPs at liquid– liquid interfaces has a crucial role in developing static and dynamic flow models that help in drug delivery and understanding the biological and physical phenomena inside the cells of the body [8].

**Citation:** Daher, A.; Ammar, A.; Hijazi, A.; Benyahia, L. Effect of Shear Flow on Nanoparticles Migration near Liquid Interfaces. *Entropy* **2021**, *23*, 1143. https://doi.org/10.3390/ e23091143

Academic Editor: Mikhail Sheremet

Received: 28 July 2021 Accepted: 25 August 2021 Published: 31 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

In addition, the dynamics of NPs in non-aqueous media, such as ionic liquids (ILs), was reported in many studies. It was shown that the design and preparation of the nanomaterials are well planned and executed, using ILs to produce tunable and functional fluid ILs-based nanomaterials, and it also was reported that ILs help to synthesize nanomaterials with various functionalized surfaces [9].

Furthermore, the non-extensivity of entropy was investigated for different sizes of colloidal Ag nanoparticles (NPs), and it was shown that the subextensivity of entropy occurs for colloidal Ag NPs. In the small size of colloidal Ag NPs and at low temperature, nonextensivity is important [10]. Taherkhani et al. used classical molecular dynamics (MD) simulations to investigate the radial distribution, glass transition, ionic transfer number, and electrical conductivity of the ionic liquid 1-ethyl-3-methylimidazolium hexafluorophosphate [EMIM][PF6] ionic liquid encapsulated in carbon nanotube (CNT), and they also studied the effect of nitrogen as a doping element in CNT on these properties of [EMIM][PF6] by MD simulation. It was shown that in the presence of nitrogen, ion transfer uses a hydrogen bonding mechanism, while in its absence, ion transfer uses a diffusion mechanism in which the cation has a significant effect on the ion transfer and electrical conductivity [11].

The behavior of NPs is strongly affected by the surrounding environmental factors and thus, external effects will modify their dynamic properties. Many researchers have studied the effects of external fields on the aggregation of NPs. The effect of shear on nanoparticle dispersion in polymer melts was investigated by Karla et al. [12], and it was shown that shear significantly slows down the aggregation of NPs and such an effect is strongly dependent on the polymer chain length and shear rate. In addition, Karla et al. [13] studied how the NPs disperse in a block copolymer system under shear flow; they found that shear can have a pronounced effect on the location of NPs in block copolymers and that it can be used as a parameter to control nano-composite self-assembly. In addition, Minh D Vo et al. [14] used dissipative particle dynamics (DPD) methods to study the effect of shear and particle shape on the physical adsorption of a polymer on carbon nanoparticles; they found that there are three possible states of the polymer adsorption on carbon nanoparticles (adsorbed, shear affected, and separated states) depending on the value of the shear rate. Besides that, the effects of shear stress on the intracellular uptake of NPs in a biomimetic microfluidic system were investigated by Kang et al. [15], and they showed that for the case of cationic NPs, as the magnitude of the shear stress increases, the intracellular uptake of such NPs maximizes at a certain value of shear stress and then decreases gradually, which ensures that the shear stress has a crucial role in various nanoparticles and drug delivery systems.

Plater et al. investigated experimentally the effect of viscosity ration of polypropylene/ Poly-e caprolactone blends on the localization of carbon black aggregates. The authors reported that the particles were dragged to the viscous phase, even when the particles were located initially into the more fluid phase, although they preferred to locate in the latter [16].

Becu et al. succeeded in visualizing a single armored droplet with nanoparticles undergoing a shear flow in another Newtonian medium. The results showed a continuous but clear slowdown of the droplet relaxation after successive strain jumps. This effect is related to the densification of the droplet interface by NPs when deformed [17].

In the current study, we focus on how shear force affects the migration of NPs to the liquid–liquid interface, which will help to understand how NPs behave under standard industrial processing conditions. The goal of this study is to implement a compact model for the simulation of the shear effect on the migration of spherical NPs near a liquid–liquid interface. Our work is based on the phase field method (PFM) for the fluids modeling, using the diffuse interface model, and molecular dynamics (MD) for modeling the motion of the nanoparticles, through which we superimpose the discrete model of NPs (using MD) on the continuum model of fluids (using PFM), which is a new idea in numerical modeling that we discussed briefly in our previous paper [18].

The content of this paper is as follows. In Section 2, we give details of the models and methods that we used in numerical simulation. We describe the diffuse interface model and give a brief description of the numerical implementation and time discretization. In Section 3, we discuss the numerical results for the migration of nanoparticles at the liquid– liquid interface. We investigate the effect of shearing on the migration of nanoparticles at the liquid–liquid interface, and we compare the results for different shear rates and different shearing times. Finally, a conclusion follows in Section 4.

#### **2. Materials and Methods**

#### *2.1. Particle–Particle and Fluid–Particle Interactions*

In this section, the discrete dynamic of particles is described. Molecular dynamics is the method of choice when one wants to study the dynamical properties of a system in full atomic detail, provided that the properties are observable within the time scale accessible to simulations. Time scale is one of the two main limitations of the method as will be discussed later. Molecular dynamics simulations are also useful when the system cannot be studied by the experimental methods mentioned above, for example, when the protein cannot be crystallized or is too big or insoluble to be studied by NMR.

To calculate the dynamics of the system, that is, the position of each atom as a function of time, Newton's classical equation of motion is solved iteratively for each atom.

Each NP is considered a rigid spherical body whose velocity (*vi*) and position (*xi*) are updated by using Newton's equation of motion, which relates the applied forces with the particle's acceleration (*ai*) according to the following equation:

$$F\_i(t) = m\_i \frac{d^2 \mathbf{x}\_i}{dt^2} = m\_i a\_i(t) \tag{1}$$

where *Fi* is the applied force on the particle *i* , and *mi* is the mass of the particle. The applied forces can be classified into particle–particle interaction forces and external forces due to the fluid (in our case, we consider drag forces, Brownian forces and shear effects).

The force on each atom is the negative of the derivative of the potential energy with respect to the position of the atom:

$$f\_{j-i} = -\nabla V(r) \tag{2}$$

If the potential energy of the system is known, then, given the coordinates of a starting structure and a set of velocities, the force acting on each atom can be calculated and a new set of coordinates generated, from which new forces can be calculated. Repetition of the procedure will generate a trajectory corresponding to the evolution of the system in time. The accuracy of the simulations is directly related to the potential energy function used to describe the interactions between particles. In molecular dynamics, a classical potential energy function is used that is defined as a function of the coordinates of each of the atoms. In macroscopic systems, the fraction of the particles near the wall is negligible, whereas in the MD simulations, this fraction is more significant, and the surface effect is of great importance. In order to reduce the surface effect and conserve the number of particles in the simulation box, periodic boundary conditions must be used so that when a particle enters or leaves the simulation region, an image particle leaves or enters this region.

The position and velocity of the NPs (*xi and vi* respectively) at each new time are calculated, using the velocity Verlet algorithm time integration scheme:

$$\mathbf{x}(t+\delta) = \mathbf{x}(t) + \mathbf{v}(t)\delta + \frac{1}{2}\mathbf{a}(t)\delta^2\tag{3}$$

$$
\upsilon(t+\delta) = \upsilon(t) + \frac{1}{2} \left[ a(t) + a(t+\delta) \right] \delta^2 \tag{4}
$$

where δ indicates the time increments for the molecular dynamics (MD).

In our model, the particle–particle interaction is calculated, using the truncated Lennard–Jones (*L Jtrunc*) potential; this potential is used in many numerical studies in order to model particle–particle interactions. See [19] for example.

$$V\_{ij\_{\rm{trunc}}}^{LJ}(r) = \begin{cases} V\_{Lf}(r\_{ij}) - V\_{Lf}(r\_c) & \text{for } r < r\_c\\ 0 & \text{for } r > r\_c \end{cases} \tag{5}$$

With

$$V\_{Lf}(r) = 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^{6} \right] \tag{6}$$

where *rc* is the cut-off distance, which is taken to be 3*σ*,  is the depth of the potential well, *σ* denotes the equilibrium distance, and *r* is center to center separation between two particles.

Thus, the interaction force acting on the *i* − *th* particle induced by the *j* − *th* particle is given by the following:

$$f\_{j \to i} = -\nabla V\_{Lf\_{\text{trunc}}}\left(r\right) \mathfrak{u}\_{j \to i} = -\frac{\partial V\_{Lf\_{\text{trunc}}}\left(r\_{ij}\right)}{\partial r\_{ij}} \mathfrak{u}\_{j \to i} \tag{7}$$

where *<sup>r</sup>ij* <sup>=</sup> *<sup>r</sup><sup>j</sup>* <sup>−</sup> *<sup>r</sup><sup>i</sup>* is the separation distance between two nanoparticles, *<sup>i</sup>* and *<sup>j</sup>* correspond to different NPs indices, and *<sup>n</sup>j*→*<sup>i</sup>* represents the unit vector that point from *xj* to *xi*.

So, for N particles, the force acting on each particle is formed by the individual interactions with all the neighboring particles:

$$F\_i = \sum\_{j \neq i}^{N} f\_{j \to i} \tag{8}$$

The external forces acting on the NPs are due the surrounding fluids and thus, each NP is affected by the Brownian force *F<sup>r</sup>* and the drag force *F<sup>d</sup>* given by the following formulas:

$$F\_r(t) = \sqrt{2 \, D \, \Delta t} \, \chi \tag{9}$$

where Δ*t* is the time step in the numerical model, *D* is the diffusion coefficient and *χ* is a normal random number whose average is zero and variance is one.

The drag effects are considered by the following:

$$F\_d = \wp \left( v\_f - v\_p \right) \tag{10}$$

where *ϕ* is the drag coefficient, *v<sup>p</sup>* is the particle's velocity and *v<sup>f</sup>* is the fluid's velocity resulting from the Navier–Stokes equation involving the concentration gradient (see [20] for details):

$$
\sigma\_f = \operatorname{div}(\nabla c \otimes \nabla \mathfrak{c} \mid \tag{11}
$$

where ⊗ corresponds to the tensor product. The definition of *c* is introduced in the next section. It denotes the phase field scalar used to describe the mixture.

#### *2.2. Diffuse Interface Model*

In the diffuse interface model, the convective Cahn–Hilliard equation is given by the following:

$$\frac{\partial \mathcal{L}}{\partial t} + \mathfrak{u}.\nabla \mathcal{c} - \nabla \mathcal{J}\_A = 0 \tag{12}$$

The diffusive flux is given by *jA* = *M*∇*μ*, where *M* denotes the mobility (scalar in the case of isotropic separation mixture, and tensor in the case of non-isotropic separation). *μ* is the chemical potential obtained from the variational derivative of the free energy with respect to the mass fraction *c*.

The Cahn–Hilliard theory [21] assumes that the driving force for diffusion is the gradient of the chemical potential, and thus, the above equation is generally written as follows:

$$\frac{\partial \mathcal{c}}{\partial t} + \mu \nabla \mathcal{c} = \nabla \cdot (M \nabla \mu) \tag{13}$$

The free energy is a double-sink function that implies that the only stable equilibrium values of *c* are +1 or −1. We use the classical symmetric form of ψ, which is used in most cases. However, other forms of non-symmetric potential can be used. The expression of free energy is assumed given by the following:

$$f(\mathbf{c}, \nabla \mathbf{c}) = -\frac{1}{2} \|\mathbf{c}\|^2 + \frac{1}{4} \|\beta \mathbf{c}^4 + \frac{1}{2} \|\nabla \mathbf{c}\|^2 \tag{14}$$

where *ε* is the gradient energy parameter and *α* and *β* are positive constants. Thus, we have the following:

$$
\mu = \frac{\delta f}{\delta c} = -ac + \beta c^3 - \varepsilon \nabla^2 c \tag{15}
$$

For this study, we apply shear forces. In order for the expression of this force to be compatible with our framework of periodic study, we consider that the shear force is taken to be periodic along the y-axis and defined by the following:

$$f\_{\text{shear}} = \frac{\sin 2 \ast \pi \ast y}{l} \tag{16}$$

This is illustrated in Figure 1.

**Figure 1.** Periodic shear force along the y-axis.

The effect of shearing is introduced into the momentum equation represented by the Navier–Stokes equation (NS) with a phase field-dependent surface force as follows [22]:

$$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = -\rho \nabla \mathbf{g} + \nabla \cdot \mathbf{n} \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T \right) + \rho \left. \mu \nabla \mathbf{c} + f\_{\text{shear}} \right| \tag{17}$$

We also consider the continuity equation for an incompressible fluid as follows:

$$
\nabla \cdot \boldsymbol{\mu} = 0 \tag{18}
$$

In the above equations, the dynamic viscosity of the fluid is denoted by η, *u* is the velocity field and *g* is the Gibbs free energy given by the following:

*g* = *f* + *p*/*ρ*, where *p* is the pressure and *ρ* is the mass density. The superscript T stands for the transpose operator.

#### *2.3. Scaling the Equations*

In order to simplify the equations and minimize the effects of round-off errors, it is preferable to use a set of dimensionless parameters. So, we scale the governing equations by defining *Uc* and *Lc* as the characteristic velocity and characteristic length, respectively, *Tc* = *Lc Uc* = *σ*. *<sup>m</sup>*  as the characteristic time and *ε<sup>c</sup>* as the characteristic energy. The characteristic length of the phase field scale is related to that of the molecular dynamics as *Lc* = 100 *σ* (fixed ratio chosen for our studies).

We introduce *r*<sup>∗</sup> = *r*/*σ*, *ε*<sup>∗</sup> = *ε*/*ε<sup>c</sup>* , *c*<sup>∗</sup> = *<sup>c</sup> cB* , *<sup>u</sup>*<sup>∗</sup> <sup>=</sup> *<sup>u</sup> Uc* , *t* <sup>∗</sup> = *tUc Lc* , *<sup>μ</sup>*<sup>∗</sup> <sup>=</sup> *<sup>μ</sup>*ξ<sup>2</sup> *<sup>ε</sup>cB* , <sup>η</sup><sup>∗</sup> <sup>=</sup> <sup>η</sup> <sup>η</sup>*<sup>c</sup>* as the normalized viscosity.

In addition, the dimensionless drag coefficient is defined as follows:

$$
\boldsymbol{\varrho}^\* = \frac{\boldsymbol{\varrho} \ \boldsymbol{\sigma}}{\sqrt{\boldsymbol{\varepsilon}\_{\boldsymbol{\varepsilon}} \ \boldsymbol{m}}}
$$

The Péclet number is defined as the product of a shear rate by a characteristic time as follows:

$$Pe = \dot{\gamma}. \ T\_c = \dot{\gamma}. \ \sigma. \ \sqrt{\frac{m}{\epsilon}}$$

ξ = *<sup>ε</sup> <sup>α</sup>* Is the interfacial thickness and *cB* <sup>=</sup> *<sup>α</sup> <sup>β</sup>* is the bulk concentration, which represents the mean field equilibrium value. Dropping the asterisk notations, we obtain the following:

$$\frac{dc}{dt} = -\mathfrak{u}.\nabla c + \frac{1}{P\_c}\,\nabla^2 \mu \tag{19}$$

$$
\mu = \mathfrak{c}^3 - \mathfrak{c} - \mathbb{C}\nabla^2 \mathfrak{c} \tag{20}
$$

$$
\nabla .\mu = 0\tag{21}
$$

$$
\nabla \mathcal{g} - \nabla . \mathfrak{v} \left( \nabla \mathfrak{u} + \nabla \mathfrak{u}^T \right) = \frac{1}{\mathbb{C} \cdot \mathbb{C}\_4} \,\mu \nabla \mathfrak{c} + f\_{\text{shear}} \tag{22}
$$

$$V\_{Lf}(\mathbf{r}) = 4\epsilon \left[ \left(\frac{1}{r}\right)^{12} - \left(\frac{1}{r}\right)^{6} \right] \tag{23}$$

From this potential, we can get the force to be the following:

$$F\_{Lf} = -\nabla V\_{Lf}\ (r)\tag{24}$$

$$F\_d = \varrho \* \left(\mathbf{v}\_f - \mathbf{v}\_p\right) \tag{25}$$

We obtain the following set of dimensionless groups: the Cahn number C and the Capillary number *Ca*, defined respectively as follows:

$$\mathbb{C} = \frac{\xi}{L\_{\emptyset}}; \,\, \mathbb{C}\_{a} = \frac{\xi \eta U\_{c}}{\rho \varepsilon \sigma\_{B}^{2}};$$

#### *2.4. Numerical Implementation*

To model the dynamics of the NPs, we use the molecular dynamics (MD) method with periodic boundary conditions implemented in order to conserve the number of NPs in the simulation box.

In addition, we use the phase field method in order to model the two fluids and the formation of the interface between them. In this method, the concentration is defined as in Equation (19).

To find the weak form, we multiply by a weighting function *c*∗ and integrate over the whole fluid domain Ω to obtain the following:

$$
\int \mathcal{c}^\* \frac{d\mathcal{c}}{dt} \, d\Omega + \int \mathcal{c}^\* \, \mu. (\nabla \mathcal{c}) d\Omega - \int \frac{1}{P\_\varepsilon} \mathcal{c}^\* \, \nabla^2 \mu \, d\Omega = 0 \tag{26}
$$

where *c*<sup>∗</sup> = *c*<sup>∗</sup> *TN* and *N* is defined as *N*<sup>1</sup> *N*<sup>2</sup> *N*<sup>3</sup> *N*<sup>4</sup> 

In the above equation *N*1, ... *N*<sup>4</sup> are the quadratic interpolation functions of the 4-node quadrilateral element:

$$\frac{dc}{dt} = \underline{N}^T \underline{\dot{c}}\tag{27}$$

.

The finite element interpolation for the gradient of the concentration is described in terms of the linear combination of the shape function derivatives, given in matrix form, by the following:

$$
\nabla \underline{c} = \underline{B} \underline{c} \tag{28}
$$

$$\mathbf{B} = \begin{bmatrix} \ N\_{1,x} & \ N\_{2,x} & \ N\_{3,x} & \ N\_{4,x} \\ \ N\_{1,y} & \ N\_{2,y} & \ N\_{3,y} & \ N\_{4,y} \end{bmatrix}.$$

Solving Equation (23), we get:

$$\underline{\mathfrak{C}}^{\*}\int \underline{\mathfrak{N}} \, \underline{\mathfrak{N}}^{T} \, d\Omega \, \dot{\underline{c}} + \underline{\mathfrak{C}}^{\*} \, \prescript{T}{}{\int \underline{\mathfrak{N}}} \, \underline{\mathfrak{u}} \, \underline{\mathfrak{B}} \, d\Omega \, \underline{\mathfrak{c}} - \underline{\mathfrak{C}}^{\*} \, \prescript{T}{}{\underline{P}\_{\varepsilon}} \int \, \mathcal{N} \, \nabla^{2} \mu \, d\Omega = 0 \tag{29}$$

This integration allows obtaining a linear system that has to be solved at each time step, which can be solved using a semi-implicit or explicit time integration scheme:

$$M\underline{\dot{\mathfrak{e}}} + G\underline{\mathfrak{e}} + F(\underline{\mathfrak{e}}) = 0 \tag{30}$$

Similarly, in order to solve the velocity equation, we can write the following:

$$
\int \mu^\* \nabla^2 \mu \, d\Omega + \int \mu^\* A \mu \nabla c d\Omega = 0 \tag{31}
$$

$$\text{where } \mathfrak{u}^\* = \underline{\mathfrak{u}}^\* \, ^T \underline{\mathcal{N}} \tag{32}$$

$$\text{and } \nabla^2 \underline{u} = \mathbb{K} \underline{u} \tag{33}$$

Solving Equation (28), we obtain the following:

$$
\underline{\mu}^\* \prescript{T}{}{\int} N \,\mathsf{K} \, d\Omega \, \underline{\mu} + \underline{\mu}^\* \prescript{T}{}{\int} N \,\mathsf{A} \, \mu \, \mathsf{B} \, d\Omega \, \underline{\mathfrak{C}} = 0 \tag{34}
$$

This integration gives the following linear system to be solved at each time step:

$$T\,\,\overline{\mu} + H\,\,\overline{\epsilon} = 0\tag{35}$$

The position and velocity of the NPs (*xi and vi* respectively) at each new time are updated, using the velocity Verlet integration scheme.

#### **3. Results and Discussion**

#### *3.1. Fluids Separation and Interface Formation*

The system is composed of two liquids that are normally immiscible. Due to some external effects (temperature for example), these two liquids may be mixed to form one thermodynamic phase. So, we can start the study from an initial state where the two fluids are totally mixed as shown in Figure 2 (left figure).

**Figure 2.** Mixture of two immiscible liquids, from the homogenous state (**left**) to the equilibrium state after the phase separation (**right**).

These figures represent a 2D plot of the simulation box representing the two fluids. The side bar in Figure 2 represents the evolution of the concentration between the two fluids in order to form the interface. The concentration is taken to vary from −1 in the first fluid (blue) to 1 in the second one (red). In our model, the area fractions of the two phases are taken, as shown in Table 1.

**Table 1.** Variation of concentrations through the medium of the two fluids.


As time goes on, and since the two liquids are normally immiscible, the molecules of each fluid immediately start to cluster together into microscopic clusters throughout the liquid. These clusters then rapidly grow and coalesce until we obtain an equilibrium state in which the two fluids are totally separated, and the interface is formed between them, as shown in Figure 2 (right).

The blue medium represents the first liquid, and the red represents the second liquid;the concentration is varying, according to the diffuse interface model. Note that our simulation box is bounded between 0 and 1 along the two axes as shown by the limiting lines in the figure, but we represent a periodic repetition of this box in Figure 2 (right) for clarity; this is done in all 2D figures throughout the paper.

#### *3.2. Introduction of Nanoparticles after the Separation of the Two Fluids and the Formation of the Interface*

After the two fluids are separated and the interface is formed between them, we introduce N nanoparticles randomly into the system.

We consider two cases regarding the concentration of NPs in the medium (total area of the NPs relative to the area of the medium of two fluids). The first case is 0.06 (200 NPs), and the second is 0.3 (1000 NPs).

#### 3.2.1. Low Concentration of Nanoparticles

Consider 200 NPs (concentration of NPs 0.06) distributed randomly within the two fluids as shown in Figure 3.

**Figure 3.** Random distribution of 200 NPs within the two fluids after phase separation.

• Neglect particle–particle interactions

For the first moment, let us neglect the interaction between the NPs via LJ potential and study the migration of NPs toward the interface. In this part, we examine the behavior of the NPs once they are introduced into the mixture of the two fluids that is already separated and the interface is formed between them. The main goal here is to study whether the external shear effect can improve the migration of NPs to the interface or not. Different cases regarding the shear effect are simulated. We start with the case involving no shear, and then we increase the Péclet number progressively; in each case, the percentage of NPs migrating to the interface is plotted.

1. No shear case.

Once introduced into the medium of the two fluids, the NPs near the interface are affected by the drag force given in Equation (7), resulting from the concentration gradient between the two fluids at the interface. These NPs are driven to the interfacial region, whereas the NPs far from the interface are not affected by this force, and thus, only 37% of the NPs migrate to the interface as shown in Figure 4.

**Figure 4.** NPs near the interface migrate to the interfacial region (**left**); percentage of NPs belonging to the interfacial region in the absence of shear (**right**).

In the simulation, time is defined as a dimensionless parameter (normalized time) t <sup>∗</sup> = tUc Lc as shown in Section 2.3 (scaling the equations).

The percentage of NPs belonging to the interface is calculated by finding the number of NPs belonging to the interfacial region relative to the total number of NPs introduced into the medium.

In this work, the interface is defined as the region of high concentration gradient, and we are able to track the position of the NPs and determine those that reach the interfacial thickness by calculating the concentration gradient at every position. By this way, we are able to determine whether each NP reaches the interface or not. It is important to note that the noise in Figure 4 (right) and all the coming figures are mainly due to the effect of the Brownian force introduced in the system. Although the effect of the other forces dominates the effect of the Brownian force, there are still some effects as seen by the percentage of error caused due to the Brownian force After the steady state is reached after introducing NPs, a shear is applied for a certain duration (T-shear), defined relative to the characteristic time. Thus, under the effect of the drag force, NPs near the interface are adsorbed. We consider different cases for which we quantify the evolution of the NPs percentage belonging to the interfacial region.

2. Simulation with shear rate = 0.4; *Pe* = 0.008, T-shear = 0.3.

In Figure 5, MD corresponds to molecular dynamics simulation, and it represents the separation of the two fluids in the absence of shear (before applying shear). This figure shows that the percentage of NPs migrating to the interface in the absence of shear is about 37%, and this percentage is not significantly affected in the case of low shear (low Péclet number *Pe* = 0.008).

**Figure 5.** Shape deformation of the two fluids under the effect of shear (**left**); percentage of NPs belonging to the interfacial region in the regions, before applying shear (green), during applying shear (red), and after stopping shearing (black), *Pe* = 0.008 (**right**).

For this case, we find that a low Péclet number does not have a noticed effect on the percentage of NPs belonging to the interface.

3. Simulation with shear rate = 7.7; *Pe* = 0.154, T-shear = 0.3.

Increasing the shear rate and thus, increasing the Péclet number to 0.154, the percentage of NPs belonging to the interfacial region increases under the effect of shear, from 37% to about 60% as shown in Figure 6.

**Figure 6.** Percentage of NPs belonging to the interfacial region in the regions, before applying shear (green), during applying shear (red) and after stopping shearing (black), *Pe* = 0.154.

4. Simulation with shear rate =15.8; *Pe* = 0.316, T-shear = 0.3.

In Figure 7, it is clear that increasing the shear rate and thus, increasing the Péclet number, enhances the migration of the NPs toward the interface to reach about 62% at the end of the shear; also, this percentage increases a little bit to reach 70% after stopping shearing.

**Figure 7.** Percentage of NPs belonging to the interfacial region in the regions, before applying shear (green), during applying shear (red) and after stopping shearing (black), *Pe* = 0.316.

This is mainly due to the fact that as the shear rate increases, more interfaces are formed in the medium (i.e., the length of the interface is increasing) and thus, the percentage of NPs belonging to the interfacial region increases. After stopping shearing, the interfaces tend to reach an equilibrium state, and they reach the separation phase. In this case, the NPs that are still near the interface are attracted to the interfacial region, due to the concentration gradient so that the percentage increases a little bit to reach 70%.

• Including particle–particle interactions.

In this part, we consider the same situation discussed above (concentration of NPs 0.06), but this time we take into account the particle–particle interactions.

Once introduced, the NPs near the interface are driven by the hydrodynamic drag force to the interfacial region, and all the NPs that are close to each other interact by the particle–particle interaction force, so the NPs form clusters around each other as shown in Figure 8.

**Figure 8.** (**a**) Random distribution of the NPs within the two fluids; (**b**) formation of clusters under the effect of particle– particle interactions.

In this part, we also study the effect of different shear rates.

1. No shear case.

Evaluating the percentage of NPs that are within the interfacial thickness, before applying any shear on the system, shows that this percentage is about 26% as shown in Figure 9. This value is smaller than that found in the case of neglecting the Lennard–Jones potential. This is mainly because the formation of clusters prevents the clustering NPs from reaching the interface. In reality, they agglomerate around the ones that are in the interfacial thickness so that the cluster is attached to the interface by some of the NPs forming it.

**Figure 9.** NPs near the interface migrate to the interfacial region (**left**); percentage of NPs belonging to the interfacial region in the regions, before applying shear (**right**).

The distribution of the NPs at the interface and the accumulation of the others over them are shown in Figure 8 above. In this case, the NPs that are attached to those at the interface are seen in the figure to belong to the first fluid (red) or the blue one (blue) and thus, they are not considered to belong to the interface. They are blocked by the ones that reached the interface before.

Now, we apply shear after the NPs near the interface are adsorbed, and we investigate the effect of shear on the migration of NPs to the interfacial region.

#### 2. Simulation with *Pe* = 0.008; T-shear = 0.3.

As for the case without the Lennard–Jones potential, small shear rates do not have an important influence on the system. Thus, the percentage of NPs belonging to the interfacial region is not significantly modified, compared to the case of no shear as shown in Figure 10.

**Figure 10.** NPs near the interface migrate to the interfacial region (**left**); percentage of NPs belonging to the interfacial region in the regions, before applying shear (green), during applying shear (red) and after stopping shearing (black), Pe = 0.008 (**right**).

3. Simulation with shear rate =15.8; *Pe* = 0.316; T-shear = 0.3.

As the shear rate increases, more interfaces are formed in the medium (i.e., the length of the interface is increasing) and thus, the percentage of NPs belonging to the interfacial region increases, as shown in Figure 11.

**Figure 11.** Percentage of NPs belonging to the interfacial region in the regions, before applying shear (green), during applying shear (red) and after stopping shearing (black), Pe = 0.316.

One can notice that the percentage in these cases differs from that in the cases discussed above, where we neglected the Lennard–Jones potential. This is mainly due to the fact that in the case when the particle–particle interactions are taken into account, the NPs accumulate around each other and form clusters at the interface.

It is noticed in the literature that embedding particles at the liquid interfaces may, for example, lead to increased stability of biphasic systems, such as Pickering's emulsions [23], or lead to a double percolation morphology followed by electrical conductivity as with

immiscible polymer blends [24]. For the latter case, several studies have focused on the competition or synergy between thermodynamics and hydrodynamics that is inerrant to mixing processes [25–28] on the particle localization. However, the observation of the particle adsorption dynamics remains rarely considered. For example, Keal et al. were able to demonstrate, by confocal microscopy tracking, the adsorption by natural diffusion of colloidal particles at a liquid–liquid interface of very low interfacial tension. To date, we are not aware of any experimental work describing the adsorption dynamics of nanoparticles at the liquid–liquid interfaces under flow.

The blue boxes (marked area) in Figure 12 show that the two interfaces come so close, and the particles accumulate between them.

**Figure 12.** Some clusters are attached between two adjacent interfaces, *Pe* = 0.316.

3.2.2. Simulation with High Concentration of NPs

As for the previous concentration, we take into account two cases. First, we neglect the Lennard–Jones (L.J.) potential and then we include it in the second case.

• Simulation of neglecting the L.J.potential.

For the initial state, consider 1000 NPs (concentration of NPs 0.3) randomly distributed within the two fluids, as shown in Figure 13a.

**Figure 13.** Random distribution of the NPs within the two fluids (**left**); more NPs migrate to the interface under the effect of shear, *Pe* = 0.316, T-shear = 0.3 (**right**).

First, neglecting the effect of the L.J. potential, we discuss the effect of different shear rates applied for the same duration (T-shear), and in addition, we study the effect of different shearing durations for the same shear rates. i.e., for the same Péclet numbers.

By increasing the Péclet number to 0.316, we find that about 66% of the NPs are driven to the interfacial thickness for shear duration (T-shear = 0.3). Thus, we again ensure that shearing enhances the migration of the NPs toward the liquid–liquid interface.

• Simulation including the L.J. potential.

In this part, we include the effect of the Lennard–Jones potential, for the case of 1000 NPs within the simulation box. We fix *Pe* = 0.3160 and T-shear = 0.3, and we examine different cases through which we vary certain parameters related to the L.J. potential.

$$1. \qquad \epsilon = 1, \sigma = 0.02$$

In Figure 14, it is clear that by increasing the concentration of NPs within the simulation box, each NP has an interaction with all the NPs close to it. Due to the high concentration of NPs, all of them are connected with each other and thus, the L.J. force is strong enough to prevent the shear from affecting the motion and the location of the NPs.

**Figure 14.** Particle–particle interactions are stronger than the effect of shear.

2. Simulation with = 0.1, *σ* = 0.02

In experimental studies, it is possible to control and vary the value of the potential between each NP and its neighbors and thus, it is possible to modify the particle–particle interactions in order to match some industrial needs [29,30].

In order to study the effect of shearing in such cases, we minimize the value of the potential well's depth, , in the definition of the L.J. potential.

The results are shown in Figure 15, where we find that in such cases, shearing affects the motion of the NPs and let them migrate toward the interfacial region. In addition, there are some NPs connected between two adjacent interfaces, and this matches the experimental results reported by Becu and Benyahia [14].

**Figure 15.** Migration of the NPs clusters to the interface under the effect of shear.

It is clear from the obtained results that by controlling the value of the potential between neighboring NPs and under the effect of shear, we are able to let the NPs form a layer on the interface, even when the concentration of the NPs in the fluids is high.

#### *3.3. Introduce the NPs Randomly into the Mixture of the Two Fluids before Phase Separation*

In this part, we study the case of NPs introduced into the medium of the homogenous mixture of the two fluids before phase separation and thus, before the formation of the interface as shown in Figure 16.

**Figure 16.** Initial state of NPs randomly distributed within the mixed fluids.

3. Simulation with no shear.

As time goes on, and since the two fluids are immiscible, they will start to separate, but this time the NPs are within the fluids and thus, their motion is affected by the phase separation.

In the absence of shear, it is clear from Figure 17 that almost all the NPs (100%) are driven to the interfacial region after a time of 2.2 relative to the characteristic time. The percentage increases progressively during the phase separation so that the NPs are affected by the force resulting from the concentration gradient throughout all the phase separation time.

**Figure 17.** Percentage of NPs belonging to the interfacial region, in the absence of shear.

Comparing these results to the case when the NPs are introduced after phase separation has occurred, it is clear that the time of NPs introduction is important.

Since the percentage of NPs adsorbing at the interface is less when the phase separation is already occurred, this seems to indicate that the potential of adsorption is higher in this case. However, when we start from the homogenous state, the potential barrier is less, and thus, more particles can be adsorbed at the interface.

In order to quantify the effect of shear on the migration of the NPs to the interface in this case, we present the results obtained for different shear rates applied to the system just at the beginning of the simulation before the phase separation starts.

4. Simulation with shear rate = 0.4; *Pe* = 0.008

As for the previous cases, a small shear rate does not have an important effect on the obtained results, compared to the case of no shear as shown in Figure 18.

**Figure 18.** Percentage of NPs belonging to the interfacial region, in the presence of small shear.


Increasing the shear rate will have a negative result with respect to the migration of NPs towards the interface. As seen in Figures 19 and 20, the percentage of NPs belonging to the interface is smaller compared to the case of no shear for all time durations.

**Figure 19.** Percentage of NPs belonging to the interfacial region in the presence of shear, *Pe* = 0.154.

**Figure 20.** Percentage of NPs belonging to the interfacial region in the presence of shear, *Pe* = 0.316.

So as seen in the figures above, we find that in the case when we introduce the NPs into the medium of the two mixed fluids before phase separation takes place, more NPs will migrate to the interface in the absence of shear.

#### **4. Conclusions**

In this study, we have investigated the effect of shear force on the migration of nanoparticles toward the interface of two immiscible liquids.

We have implemented a numerical simulation through which we modeled the two fluids using the diffuse interface model. Two cases were investigated. The first one is to introduce the NPs randomly into the medium of two fluids that have been separated and the interface was formed between them. The second case is to introduce the NPs randomly into the mixture of the two fluids before phase separation takes place. For the first case, if we leave the fluids without any external effect, we find that only a small percentage, not greater than 30% of the NPs will migrate towards the interface. These NPs are the ones that are close to the interface when we introduce them randomly into the medium. Their migration toward the interface is mainly due to the effect of the drag force related to the concentration gradient in the interfacial region.

We have found that inducing a shear onto the system enhances the percentage of NPs belonging to the interface. It has been shown that there are two factors affecting the percentage of NPs migrating towards the interface, the value of the shear rate modified using the Péclet number, *Pe and* the shear duration T-shear. Introducing a shear defined by *Pe* = 0.154 for a duration of 0.63 drives 66% of the NPs to the interface compared to 30 % for the case of no shear effect. The same result was obtained using *Pe* = 0.316 for duration of 0.3. This ensures that both factors play a crucial role regarding the migration of NPs towards the interface.

In addition, we have discussed the effect of including or neglecting the particleparticle interaction using the Lennard–Jones potential and we have found that for low concentration of NPs, some differences appear regarding the percentage of NPs within the interface. This is mainly because due to particle–particle interactions NPs close to each other will form clusters and accumulate around each other, so that then NPs first attached to the interface will attract the close ones to accumulate around them so that we find clusters of NPs attached to the interface. In addition, for the same case we find that some clusters will be attached between two close interfaces which match with some experimental observations [14]. On the other hand, we have found that for the case of high concentration of NPs there are several observations depending on the value of the L.J. potential well. For  = 1 the particle–particle interaction will be strong enough and will compete the effect of shear and thus the NPs will not be driven to the interface. However, for  = 0.1 the NPs will be driven to the interface, where they form a layer at the interfacial region, with some of them are connected with other NPs on adjacent interfaces. From the experimental point of view, it is well known that we can modify the value of the particle–particle potential interactions in order to match some industrial needs and thus the NPs can be driven under the effect of shear to the interface even though when their concentration in the fluids is high. So for this case, we have shown that shearing is a key factor in enhancing the migration of the NPs towards the interface.

On the other hand, we have found that NPs introduced randomly into the mixture of the two fluids before phase separation takes place will migrate faster and in a higher percentage to the interface in the absence of shear.

It is clear that the time of NPs introduction is important. Since the percentage of NPs adsorbing at the interface is less when the phase separation already occurred, seems to indicate that the potential of adsorption is higher in the case. Whereas, when we start from the homogenous state, the potential barrier is less, thus, more particles can be adsorbed at the interface.

These results help us to understand how NPs behave under standard industrial processing conditions. We obtain more information regarding new methods for synthesizing nanomaterial films; in addition, these results also help to understand the behavior of small NPs within the cells of organisms where they are greatly affected by the flow rate of the surrounding fluids.

**Author Contributions:** Data curation, A.D.; Formal analysis, A.D., A.H. and L.B.; Funding acquisition, A.H.; Investigation, A.D.; Methodology, A.D.; Project administration, A.H., A.A., Software, A.D. and A.A.; Supervision, A.A. and A.H.; Validation, A.D.; Writing—original draft, A.D.; Writing review & editing, A.D. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by grant (2020) from the Lebanese University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

