Shear-Induced Migration of Rigid Particles near an Interface between a Newtonian and a Viscoelastic Fluid.

Simulations of rigid particles suspended in two-phase shear flow are presented, where one of the suspending fluids is viscoelastic, whereas the other is Newtonian. The Cahn-Hilliard diffuse-interface model is employed for the fluid-fluid interface, which can naturally describe the interactions between the particle and the interface (e.g., particle adsorption). It is shown that particles can migrate across streamlines of the base flow, which is due to the surface tension of the fluid-fluid interface and a difference in normal stresses between the two fluids. The particle is initially located in the viscoelastic fluid, and its migration is investigated in terms of the Weissenberg number Wi (shear rate times relaxation time) and capillary number Ca (viscous stress over capillary stress). Four regimes of particle migration are observed, which can roughly be described by migration away from the interface for Wi = 0, halted migration toward the interface for low Wi and low Ca, particle adsorption at the interface for high Wi and low Ca, and penetration into the Newtonian fluid for high Wi and high Ca. It is found that the angular velocity of the particle plays a large role in determining the final location of the particle, especially for high Wi. From morphology plots, it is deduced that the different dynamics can be described well by considering a balance in the first-normal stress difference and Laplace pressure. However, it is shown that other parameters, such as the equilibrium contact angle and diffusion of the fluid, are also important in determining the final location of the particle.

Introduction

Understanding the dynamics of rigid particles in viscoelastic two-phase flows is of great importance for many technological applications. For example, in the field of polymer processing, immiscible polymers are combined to create novel materials, to which rigid particles are often added to improve toughness.[1] During the processing of these materials, distinct domains are formed by the immiscible molten polymers, with a fluid–fluid interface between them. The final material properties will depend to a large degree on the location of the particles in the material (e.g., in one of the polymers and/or at the interface). Polymer blends with rigid fillers have therefore received considerable attention over the last few decades.[2−5] If particles are located at the interface between the fluids, then stable emulsions can be formed, known as Pickering emulsions.[6] Particles at interfaces can self-assemble in interesting patterns, strongly influencing the rheology of the fluid–fluid interfaces.[7] Moreover, it was shown that the coalescence of droplets in polymer blends can be significantly delayed by particles at the fluid–fluid interface.[8] Due to the high viscosity of polymeric fluids, the Brownian motion of the particles is, in general, negligible; therefore, flow is essential to moving particles close enough to the interface such that they are adsorbed at the interface. Under the influence of flow, particles can migrate toward one of the polymeric domains, as was observed by Elias et al. for non-Brownian aggregates of nanosilica particles in polymer blends under shear.[9] One of the parameters that influences the location of the particles is the surface chemistry of the particles, which can favor one phase more than the other. However, the particles still need to come close enough to the interface to “feel” the other phase, thus in the absence of Brownian motion, surface chemistry alone cannot explain the migration. Elias et al. explain the migration of particles by collisions between particles and droplets due to flow. In this article, we investigate an alternative explanation of the migration of particles in two-phase viscoelastic flows, which is based on a difference in normal stresses in the suspending fluids.

When inertia is neglected, a single rigid spherical particle suspended in a Newtonian fluid under shear will remain on the streamline of the base flow. However, particles can migrate across streamlines in Newtonian shear flows if inertia does play a role.[10] Furthermore, both experiments[11−13] and simulations[13−15] show that particles can migrate across streamlines if the suspending fluid is viscoelastic, even if inertia is negligible. The migration of particles in creeping viscoelastic flows can be attributed to gradients of normal stresses. For example, in a wide-gap Couette flow device, the shear rate near the inner cylinder is higher than that near the outer cylinder, which leads to a gradient in normal stresses, moving the particle toward the outer cylinder.[13] Gradients of normal stresses also occur if the material properties vary across the domain, even if the shear rate is constant. This idea motivated the current study which entails the numerical investigation of rigid particles suspended in two-phase fluids, where one of the fluids is Newtonian and the other is viscoelastic. The particles are initially located in the viscoelastic fluid, and the effect of the rheology of the fluids and the surface tension on the migration of the rigid spheres is investigated.

Theoretical Model

The problem considered in this article is a rigid spherical particle of radius a in a two-phase shear flow, as depicted in Figure . The upper fluid is Newtonian, whereas the lower fluid is viscoelastic. Between the fluids there exists a fluid–fluid interface which is endowed with a surface tension σ. Effects from gravity are neglected (i.e., it is assumed that the particle and fluids have similar densities). The top and bottom walls move in the x direction with respective velocities of Uw and −Uw, yielding a global shear rate of γ̇ = 2Uw/H. Initially, the particle is located in the lower fluid, and the motion in the y direction due to the imposed shear flow is investigated.

Figure 1

Problem consisting of a spherical particle in a two-phase shear flow. The fluid–fluid interface is shown as a blue surface, whereas the rigid walls are depicted as gray surfaces. The domain is bounded in the y direction by rigid walls, and periodicity is assumed in the x direction. For the z direction, symmetry is assumed in the xy planes at z = 0 and z = W/2. (The origin is located at the centroid of the rectangular box.)

Cahn–Hilliard Theory

To describe the fluid–fluid interface, a diffuse-interface model is employed that can naturally describe singular events such as droplet coalescence[16] and moving contact lines.[17] A phase-field variable ϕ is introduced, which attains constant values inside each fluid but varies continuously across the interface. The variable ϕ can thus be identified with the local composition of the fluid. The evolution of ϕ is described by the Cahn–Hilliard equation[18]where D()/Dt is the material derivative, M is the mobility, which is assumed to be constant in this article, and μ is the chemical potential. To obtain an expression for the chemical potential in terms of the composition, the total free energy of the system is written as an integral over the volume Ω and physical boundaries Γwhere F is the total free energy, f is the free energy density, and fw is the wall free energy. The assumption of Cahn and Hilliard is that the free energy f depends on the local composition ϕ and on gradients of ϕ, adding weak nonlocal interactions[18]where the first term on the right-hand side represents a double-well potential that can be scaled by χ and which has minima at ϕ = ±1 (i.e., the bulk values of the phase-field parameter). Moreover, κ is the gradient energy parameter. The second term in eq is minimized when the gradients of ϕ vanish, thus this term promotes mixing. The combination of the two terms in the expression for the free energy leads to a diffuse interface. An expression for the chemical potential in the bulk can be obtained by taking the variational derivative of the free energy with respect to the compositionwhere δ()/δϕ is the variational derivative. For a planar interface in equilibrium (ϕ = ϕ(x)), eqs and 4 can be solved analytically to yield the interface profile ϕ(x) = tanh[x/(21/2ξ)], where is a measure of the interface thickness. The inclusion of gradients of ϕ in the free-energy expression leads to the excess energywhere σ* is the Cahn–Hilliard interfacial energy, which will be identified with the interfacial energy σ. However, it should be emphasized that flow can move the interface profile away from its equilibrium solution, yielding local variations in the interfacial energy.[19]

The wall free energy is a function of the composition at the wall and can be used to include fluid–solid interfacial tensions.[20] The wall free energy is given bywhere ζ is the (scalar) wetting potential. The function fw is monotonically increasing (for ζ > 0) or decreasing (for ζ < 0) for −1 ≤ ϕ ≤ 1 and can thus be used to impose a different fluid–solid interfacial tension for ϕ = 1 than for ϕ = −1. This difference in the fluid–solid interfacial tension can be used together with Young’s equation to yield a relation between the wetting potential and the equilibrium contact angle θc: ζ = 4σ* cos θc/3. A boundary condition for ϕ can be obtained by considering variations in fw at the boundary, which leads to[21]where ∂()/∂n is the directional derivative in the direction normal to the boundary. Note that it is possible to extend this model to include a nonequilibrium (dynamic) contact angle.[22] However, for simplicity in the analysis, only the boundary condition as given in eq is used, which imposes the equilibrium (static) contact angle on the physical boundaries Γ. To obtain a boundary condition for the chemical potential, a no-flux condition on all physical boundaries is assumed: ∂μ/∂n = 0. To integrate the Cahn–Hilliard equation in time, an initial condition is needed for the composition field ϕ. For this, we use the equilibrium tanh interface profile, with the interface location (defined by ϕ = 0) as shown in Figure .

Flow Equations

Since our main interest is in highly viscous polymeric fluids, it is assumed that inertia does not play a role (i.e., creeping flow is considered). Furthermore, the fluids are incompressible and density-matched, which yields the balance of momentum and balance of masswhere σ is the Cauchy stress tensor and is the fluid velocity. Due to the assumption of creeping flow, the particles are force- and torque-free, which is expressed aswhere ∂P is the particle boundary, is the outwardly directed unit normal on the particle boundary, and = [X, Y, 0] is the location of the center point of the particle (where the z coordinate remains zero due to the symmetry assumption as shown in Figure ).

The Cauchy stress tensor used in this article is given bywhere p is the pressure, is the unit tensor, τs is the Newtonian stress tensor, τp is the viscoelastic (polymer) stress tensor, and τc is the capillary stress tensor (i.e., the surface tension).

To describe the two fluids, we assume that the material parameters are a function of the local composition ϕ using a linear mixing rule. The initial composition field is chosen such that ϕ = 1 in the upper (Newtonian) fluid and ϕ = −1 in the lower (viscoelastic) fluid. The Newtonian stress is given bywhere ηn is the viscosity of the Newtonian fluid and ηs is the solvent viscosity of the viscoelastic fluid. Note that a solvent in the viscoelastic fluid is mainly included for (numerical) stability and is chosen to be much smaller than ηn. For the viscoelastic fluid, the Giesekus model is employed, which captures many of the essential rheological features of polymeric fluids (i.e., shear-thinning and strain-hardening).[23] It is assumed that the modulus G is a function of ϕ such that it vanishes in the Newtonian fluid, which implies a linear decrease in polymer density with ϕ across the interface.[24] The polymer stress is written aswhere is the conformation tensor whose evolution is described bywhere λ is the relaxation time, is the upper-convected derivative given by , and α is a parameter that controls the shear-thinning behavior of the fluid. The zero-shear viscosity of the viscoelastic fluid is given by η0 = Gλ + ηs. To integrate eq in time, an initial condition is needed for the conformation tensor, for which the stress-free state is assumed: (t = 0) = .

Finally, the capillary stress is considered. In the Cahn–Hilliard framework, the capillary stress arises naturally due to the addition of the ∇ϕ term to the expression for the free energy of the fluid (see eq ). Using a variational approach, this stress tensor is found to be[24−26]where an isotropic term was added to ensure that the stress is parallel to the interface.[25,26] The capillary stress as defined in eq was found to have superior numerical convergence properties in the presence of freely floating particles, as was shown in ref (27) and will be used for all simulations presented in this article. Moreover, it can be shown that the stress tensor as written in eq converges to the sharp-interface stress tensor in the limit of a small interface thickness.[19]

On all physical boundaries, a no-slip condition is assumed for the fluid velocity. On the top and bottom walls (as shown in Figure ), a velocity is imposed in the positive and negative x directions, respectively, with a magnitude Uw. Moreover, the fluid velocity on the particle boundary satisfies = + ω × ( – ), where = [U, U, 0] is the translational velocity of the particle and ω is the angular velocity of the particle. Note that our numerical approach ensures that the particle velocities are such that the conditions given in eqs and 11 are satisfied. Finally, the motion of the particle is described by the following kinematic relationwhere the rotation angle of the particle does not need to be updated due to the particle being spherical. To integrate eq in time, the initial particle position is needed, which is placed at a distance of 2 times the particle radius below the interface, as shown in Figure .

Numerical Method

To solve the Cahn–Hilliard equation, the mass and momentum balance, and the evolution equation for the conformation tensor, we use the finite element method with adaptive meshing and adaptive time stepping. Due to the symmetry and periodicity assumptions (see Figure ), only half of the particle has been simulated. More information on the numerical method can be found in the Supporting Information or in our previous work on particles in Cahn–Hilliard fluids in Newtonian flows[27,28] and viscoelastic flows.[19]

Results

Dimensional Analysis

All results will be presented in dimensionless form, and the relevant dimensionless groups are introduced in this section. Even though all simulations performed for this article use the Cahn–Hilliard framework to describe the fluid–fluid interface, for clarity we will first treat the dimensionless groups that arise without considering the interface thickness and diffusion of the fluids. From dimensional analysis, we findwhere Wi is the Weissenberg number, Ca is the capillary number, β is the ratio between the solvent viscosity and zero-shear viscosity of the viscoelastic fluid, and δ is the ratio between the viscosity of the Newtonian fluid and the zero-shear viscosity of the viscoelastic fluid. The results will be presented for varying Wi and Ca. The viscosity ratios are set to δ = 1 and β = 0.1 for all simulations, whereas the Giesekus mobility is set to α = 0.2 for all simulations. The size of the domain will be fixed at L = W = 4a and H = 40a, and the particle is initially located in the viscoelastic fluid, with its center point at a distance of 2a from the interface (see Figure ). Note that the periodic and symmetry boundary conditions (see Figure ) imply that the system under consideration is actually an array of spheres migrating simultaneously. By using L = W = 4a, interactions between the particles are likely to play a role, but this domain size was chosen to keep the problem numerically tractable. Moreover, we believe that this problem is of high practical relevance since in polymer processing particle volume fractions are typically high and strong particle–particle interactions are to be expected. The distance in the y direction was chosen to be large enough that the walls have a minimal influence on the migration behavior of the particle, as shown in the Supporting Information. For a better interpretation of the results, the relative viscosity and the relative first-normal stress coefficient of the viscoelastic suspending fluid are presented in Figure . Furthermore, we note that the local shear rates in the Newtonian and Giesekus fluids differ from the global shear rate γ̇ due to the shear-thinning behavior of the Giesekus fluid. To aid in the analysis, the expected velocity profile and local shear rates for varying Wi are shown in Figure . These results were obtained using a sharp-interface model with a continuity of traction and velocity at the fluid–fluid interface and imposed velocity at the walls, similar to those in reference (19).

Figure 2

Relative viscosity ηr = σ/(η0γ̇) (a) and the relative first-normal stress coefficient Ψ1r = N1/(η0λγ̇2) (b) for a Giesekus fluid in steady shear with α = 0.2 and β = 0.1.

Figure 3

Expected velocity profile and local shear rate for a two-layered shear flow, where the upper fluid is Newtonian and the lower fluid is a Giesekus fluid with α = 0.2, β = 0.1, and δ = 1.

The dimensionless groups presented in eq , together with the size of the domain and the initial location of the particle and interface, completely govern the problem in case a sharp-interface approach is used. However, the use of the Cahn–Hilliard framework to describe the two fluids adds additional parameters: the interface thickness ξ, the mobility M, and the equilibrium contact angle θc, which yields three additional dimensionless groupswhere Cn is the Cahn number, S represents a ratio between the diffusion length and the macroscopic length with being an effective viscosity,[20] and θc is the equilibrium contact angle. To estimate Cn and S, an order-of-magnitude estimation is performed using physical values. In the diffuse-interface model, different definitions can be used to define the interface thickness. Here, the definition as commonly used in experiments is adopted, which defines the interfacial thickness as the distance between the intersections between the gradient of ϕ at ϕ = 0 and the bulk values of ϕ, which is given by 2(2)1/2ξ (see Figure ). For fluids consisting of small molecules, the interface thickness is on the order of 1 nm,[29] but macromolecular fluids often have interfaces that are thicker;[30,31] therefore, an interface thickness of 10 nm is assumed to estimate Cn. The particle size is set to a diameter of 100 nm, which was given as the typical size of non-Brownian aggregates of nanosilica particles in ref (9). These assumptions yield a Cahn number (as defined in eq ) of Cn = 0.071, which will be used for all simulations in this article. For the constant mobility M, typical values found in the literature are M = 10–17 m5 s–1 J–1.[29,32] Assuming a value for the viscosity of 30 Pa·s, which is typical for a viscoelastic wormlike micellar surfactant solution,[33] leads to a value of S = 0.34. However, the value of S is expected to be lower in macromolecular fluids, where diffusion typically occurs much more slowly.[34] Furthermore, measuring the Cahn–Hilliard mobility M is challenging, especially close to solid boundaries. Therefore, a value of S = 0.1 is used as a base case. Since the mobility can play a large role in phase-field simulations,[35] the influence of varying S is investigated as well. The contact angle θc is relevant only for the physical boundaries, which are the particle boundary and the top and bottom walls, as shown in Figure . On the particle boundary, the boundary condition given by eq is imposed for varying values of θc to investigate the influence of the contact angle on the migration of the particle. The fluid–fluid interface remains far from the top and bottom walls, and thus the contact angle does not play a role there.

Figure 4

Interface thickness is defined as the distance between the intersections of the gradient of ϕ at ϕ = 0 and the bulk values of ϕ, yielding a value of 2(2)1/2ξ.

Four Regimes of Particle Migration

By performing simulations for varying values of Wi and Ca, we identified four possible scenarios for particle migration:In Figure , results are presented for Ca = 1 and S = 0.1, and by changing only Wi, all of the scenarios are reproduced. For Wi = 0, the particle clearly moves downward, away from the interface. As Wi is increased, the particle moves toward the interface, but this motion is halted for Wi = 1. The particle makes contact with the interface for both Wi = 2 and Wi = 3, which is accompanied by a rapid increase in the vertical velocity of the particle. For Wi = 3, the migration velocity then decreases, and the particle remains attached to the interface. For Wi = 2, the particle keeps moving through the interface and detaches into the upper fluid. In the next subsections, we will look at each scenario in more detail.

Figure 5

Particle vertical position (a) and particle velocity (b). Ca = 1, θc = 90°, and S = 0.1.

the particle migrates away from the interface;

the particle migrates toward the interface, but its migration is halted;

the particle penetrates the interface into the Newtonian fluid; or

the particle migrates toward the interface and is adsorbed at the interface.

Migration Away from the Interface

The first regime occurs when the normal stress differences in the lower fluid are absent. Due to the deformation of the interface, a Laplace pressure will be build up which effectively pushes the particle downward, as shown in Figure . Note that the deformation of the interface is small and hardly visible yet large enough to yield a negative vertical particle velocity. The particle location Y and particle velocity U as a function of strain are presented in Figure for Wi = 0 (both fluids are Newtonian) for several values of Ca. Oscillations in the particle velocity are observed, which are due to the initial disturbance of the interface and are rapidly smoothed out by the surface tension. A negative migration velocity is observed for all values of Ca, with a magnitude that becomes larger for smaller Ca. The migration velocity decreases as the particle moves away from the interface but does not decrease to negligible values within the duration of the simulation.

Figure 6

Snapshots of the particle migrating away from the fluid–fluid interface (represented by the blue surface). Wi = 0, S = 0.1, Ca = 1, and θc = 90°.

Figure 7

Particle vertical position (a) and particle velocity (b). Wi = 0, S = 0.1, and θc = 90°.

Halted Particle Migration

As Wi is increased, a migration toward the Newtonian fluid can be induced. As shown in Figure , the flow around the upward-moving particle leads to a deformation of the interface, which in turn yields a positive Laplace pressure on the lower side of the interface. This Laplace pressure effectively pushes the particle down, and if the surface tension is large enough, the migration toward the Newtonian fluid can be halted. In Figure , the particle location and velocity are presented for Wi = 0.5 and several values of Ca. It can be observed that the particle indeed attains a positive velocity in the y direction immediately: it is migrating toward the interface. The magnitude of the migration velocity depends on the capillary number, and faster migration is observed for higher Ca. This can be readily explained by an increase in the surface tension, leading to larger Laplace pressures exerting a downward force on the particle. Around tγ̇ = 100, the migration velocity goes through a maximum, after which it decays to values that are several orders of magnitude lower than the initial migration velocity. For Ca = 0.5, the migration velocity appears to decrease to zero following a power law, whereas the other values of Ca show that the migration velocity reaches a small but finite value. In all cases, the particle reaches a stable position below the interface, but the particle gets closer to the interface as Ca is increased, as can be seen in the particle location. Again, this effect can be readily explained by considering a balance between the normal stresses and the Laplace pressure: for increasing Wi, the normal stresses are large, and a larger deformation of the interface is needed to yield the necessary Laplace pressure to halt the migration of the particle.

Figure 8

Snapshots of halted migration (the fluid–fluid interface is represented by the blue surface). Wi = 1, S = 0.1, Ca = 1, and θc = 90°.

Figure 9

Particle vertical positions (a), particle velocity (b), and particle velocity on a log scale (c). Wi = 0.5, S = 0.1, and θc = 90°.

The trace of the polymer stress is shown for Wi = 1 and Ca = 1 in Figure , where it can be clearly observed that the polymer stresses exist only in the lower fluid. Moreover, the viscoelastic stresses underneath the particle remain relatively constant as the particle is migrating. As the particle moves upward, the layer of fluid in between the particle and the interface decreases, which might explain the initial increase in migration velocity: normal stresses act in the viscoelastic fluid both above and below the particle, but due to the presence of the Newtonian fluid, the normal stresses will be higher below the particle. As the particle moves upward, the layer of viscoelastic fluid above the particle becomes thinner, further decreasing the normal stresses above the particle, yielding an increase in migration velocity. It can be observed that large areas of stress exist below the particle, which remain relatively constant as the particle is moving toward the interface.

Figure 10

Snapshots of the trace of the polymer stress for halted migration. Wi = 1, S = 0.1, Ca = 1, and θc = 90°.

Penetration through the Interface

By increasing Wi further, the normal stresses in the lower fluid can be increased to such an extent that the interfacial tension is not large enough to halt the migration of the particle. As a result, the particle makes contact with the interface and can subsequently move into the upper fluid, a process that can be simulated by the use of a diffuse-interface model for the fluid–fluid interface. For Wi = 2 and Ca = 1, the particle migrates into the upper fluid, as shown in Figure . As can be observed in Figure , a sudden and large increase in migration velocity is observed as the particle makes contact with the interface and when the particle detaches from the interface. The trace of the viscoelastic stress is shown in Figure , where it can be observed that the stresses in the lower fluid are higher compared to the case of Wi = 1 (as presented in Figure ). It can furthermore be observed that the interface moves along the particle boundary but not necessarily with the rotation of the particle. (Note that the particle is rotating in the clockwise direction.) This contact-line motion is mainly governed by S, which therefore is likely to have a large influence on the particle dynamics, especially when the particle has made contact with the interface. The influence of S on the particle migration will be further investigated in one of the following sections.

Figure 11

Snapshots of the particle penetrating the fluid–fluid interface and migrating into the Newtonian fluid (the fluid–fluid interface is represented by the blue surface). Wi = 2, S = 0.1, Ca = 1, and θc = 90°.

Figure 12

Snapshots of the trace of the polymer stress for a particle migrating into the upper fluid. Wi = 2, S = 0.1, Ca = 1, and θc = 90°.

Adsorption at the Interface

The final scenario that can occur is the adsorption of the particle at the interface. Similar to the interface penetration regime, the particle makes contact with the interface, but in this case the particle remains attached to the interface. Snapshots of this scenario are presented in Figure . The initial dynamics are similar to the penetration scenario, with the particle making contact with the interface and the interface subsequently moving along the particle boundary. However, in this case the downward force of the interface pulling on the particle is enough to balance the upward force that arises due to the gradients in normal stresses. The viscoelastic stress is shown in Figure for Wi = 3 and Ca = 1, where again larger stresses can be observed compared to Wi = 1 and 2. It is interesting that these stresses are large enough to push the particle into the Newtonian fluid for Wi = 2, whereas the particle remains attached to the interface for Wi = 3. A possible explanation can be found in the angular velocity of the particle which is known to decrease with increasing Wi.[36,37] The angular velocity around the z axis (denoted by ω and defined as positive in the clockwise direction) is shown for Ca = 1 and varying Wi in Figure , where it can be seen that the angular velocity indeed decreases with increasing Wi. As the particle is rotating while being adsorbed at the interface, the contact line of the fluid–fluid interface with the particle boundary slips by means of Cahn–Hilliard diffusion. For lower angular velocities, the contact line can remain on the particles, whereas the particle spins off the interface for higher angular velocities, explaining why the particle remains more easily attached at the interface for higher Wi. Note that an isolated particle rotates with an angular velocity of ω/γ̇ = 0.5 in a Newtonian fluid.[38] However, due to interactions between the periodic particles, the angular velocity is slightly lower for Wi = 0. Furthermore, for Wi = 2 the particle is located in the Newtonian fluid for tγ̇ > 300, but the angular velocity appears to be smaller than for Wi = 0. This can be readily explained by the shear rate being lower in the Newtonian fluid for Wi = 2, as shown in Figure . When using the local shear rate in the Newtonian fluid to scale the angular velocity, a similar value is found for Wi = 2 as for Wi = 0.

Figure 13

Snapshots of the particle being adsorbed at the fluid–fluid interface (represented by the blue surface). Wi = 3, S = 0.1, Ca = 1, and θc = 90°.

Figure 14

Snapshots of the trace of the polymer stress for a particle being adsorbed at the fluid–fluid interface. Wi = 3, S = 0.1, Ca = 1, and θc = 90°.

Figure 15

Particle angular velocity for Ca = 1, S = 0.1, and θc = 90°.

Morphology Plots

Having established the different scenarios of particle migration that can take place, we will now proceed with the analysis of the final location of the particle. This is done by evaluating the position of the particle at tγ̇ = 500 and creating morphology plots that show the final location of the particle. In Figure a, a morphology plot is presented for S = 0.1 and θc = 90° and varying Wi and Ca. Distinct regions can be observed where the different scenarios take place. These can be roughly described by migration away for Wi = 0, halted migration for low Wi and low Ca, particle adsorption for high Wi and low Ca, and penetration for high Wi and high Ca.

Figure 16

Morphology plot for S = 0.1 and θc = 90°. Each point denotes the location of the particle at tγ̇ = 500. The global Weissenberg number Wi is used in (a), whereas the local Weissenberg number Wî is used in (b).

As explained earlier, the migration behavior is attributed to differences in normal stresses between the two fluids. Using this idea, the particle location can be predicted by a simple force balance on the particle. The relevant stresses are the first-normal stress difference of the viscoelastic fluid, defined in steady simple shear by N1 = σ – σ, and the Laplace pressure due to a curved interface is given by Δp ∼ σ/a, where the curvature is assumed to be proportional to the particle radius (see Figure ). Furthermore, to describe the elastic stresses properly, we introduce the local Weissenberg number Ŵi, which is defined using the local shear rate in the viscoelastic fluid (see Figure ). The final location of the particle as a function of Ca and Wî is shown in Figure b. In the same figure, the isoline of N1/Δp = 1 is plotted, which is found to give a good description of the area where particle penetration takes place. These results support the idea that for S = 0.1 and θc = 90°, particle penetration is mainly governed by a balance between the first normal stress difference and Laplace pressure.

Influence of the Contact Angle

Next, the influence of the equilibrium contact angle on the final location of the particle is investigated. Similar to the morphology plot as presented in Figure , simulations were performed for θc = 60 and 120° (S = 0.1), and the final location was marked in a morphology plot. These morphology plots are presented in Figure a for θc = 60° and in Figure b for θc = 120°. Note that the contact angle is measured through the lower fluid, and thus a contact angle smaller than 90° means that the particle favors the lower fluid, whereas a contact angle larger than 90° means that the particle favors the upper fluid. It can clearly be observed that the contact angle plays a large role in determining the final location of the particle. For θc = 60°, the region where penetration takes place is smaller compared to that for θc = 90°. The opposite is observed for θc = 120°, where a significant increase in the penetration region can be seen. The region where particle adsorption takes place is largest for θc = 60°. These results indicate that the most likely route to getting particles at the interface by means of elastically induced migration is by placing them in the fluid which has both the highest N1 and favorable chemistry for the particle surface.

Figure 17

Morphology plot for varying contact angle (S = 0.1). Each point denotes the location of the particle at tγ̇ = 500.

Influence of the Cahn–Hilliard Mobility

We conclude by investigating the role of the mobility of the fluids. A morphology plot is shown for θc = 90° and S = 0.01 in Figure , where it can be seen that the halted migration area is increased significantly compared to that for S = 0.1. Moreover, the particle adsorption region is completely absent. The increase in the area where halted migration takes place can be explained by a compression of the tanh profile of the interface, which cannot be compensated for due to the low value of S. This compression of the interface yields a local increase in the effective surface tension,[19,24] halting the migration of the particle at an earlier stage. Moreover, a decrease in S implies a decrease in the diffusion length (as shown in its definition in eq ). As the particle makes contact with the interface, the diffusion-governed slip of the interface is much less pronounced, causing wetting failure as the particle is rotating at the interface, which explains the absence of the adsorption region. We conclude by showing snapshots of a particle penetrating the interface for lower values of S in Figure . Due to the lower diffusion of the fluids, droplets of the lower fluid can remain attached to the particle. Particles that are completely enclosed by the lower fluid, as observed in the 2D simulations in reference (19), are not observed in the 3D simulations presented in this article. Possibly, lower values of Cn are necessary for this, which is outside the scope of this article.

Figure 18

Morphology plot for S = 0.01 and θc = 90°. Each point denotes the location of the particle at tγ̇ = 500.

Figure 19

Snapshots of the particle penetrating the fluid–fluid interface and migrating into the Newtonian fluid. (The fluid–fluid interface is represented by the blue surface.) Wi = 1, S = 0.01, Ca = 4, and θc = 90°.

Discussion and Conclusions

We have presented simulations of particle migration in two-phase shear flow, where one of the fluids is viscoelastic and the other is Newtonian. The Cahn–Hilliard diffuse-interface model is used to describe the two fluids, and viscoelastic stresses are present in only one of the fluids. Initially, the particle is located in the viscoelastic fluid, but the particle has a tendency to migrate toward the Newtonian fluid due to the shear flow. We believe that this is caused by gradients in normal stresses, similar to the migration of particles toward the outer cylinder in a wide-gap Couette device.[13] However, for particles in two-phase viscoelastic/Newtonian shear flow, the gradients of normal stresses are due to inhomogeneous material parameters instead of an inhomogeneous shear rate. The results indicate that a force balance based on the first-normal stress difference of the viscoelastic fluid and the Laplace pressure can be used to predict the penetration of the particle into the Newtonian fluid. However, both the contact angle of the fluid–fluid interface with the particle boundary and the diffusion of the fluids play a large role in determining the final location of the particle. Furthermore, it was shown that the angular velocity of the particles (that is known to decrease with increasing Weissenberg number[36,37]) determines to a large extent if a particle remains adsorbed at the interface. The model can be easily adapted to simulate multiphase viscoelastic fluids by including multiple viscoelastic modes. An interesting question one could ask is whether particle migration near an interface between two viscoelastic fluids is governed by the difference in the first-normal stress difference of the two fluids. This will be a topic of future research.

The thickness of the diffuse interface was estimated using physical values, where a particle diameter of 100 nm was used, similar to the experiments presented by Elias et al.[9] Future investigations will include the influence of the interface thickness on the migration of particles near fluid–fluid interfaces (or, similarly:, changing the particle size), possibly with relation to a sharp-interface model. In the present model, the motion of the contact line across the surface of the particle is governed by the Cahn–Hilliard mobility in a phenomenological sense: contact line pinning and hopping, which are known to be important in the adsorption of particles at interfaces,[39] are not described explicitly. A possible extension of the model is to explicitly describe the roughness of the particle surface. Moreover, detailed experimental results on particles interacting with fluid–fluid interfaces in (viscoelastic) flows are crucial to verifying the model.

The numerical model used in this article was set up in a general fashion and can easily be adapted to simulate other cases of particles in (viscoelastic) multiphase flows. For example, by changing the particle shape, type of flow, and rheology of the suspending fluids, many interesting problems that are of practical relevance can be studied. Some applications that one might think of is the smart design of materials by directing particles to a certain fluid or to the fluid–fluid interface or using the rheological properties of the suspending fluids to control the motion of particles in microfluidics.[40]