Direct numerical simulations of emulsion flows

Direct numerical simulations of emulsion flows

F. R. Cunha^I; M. H. P. Almeida^II; M. Loewenberg^III

^IUniversidade de Brasília, Departamento de Engenharia Mecânica, 70910900 Brasília, DF. Brazil. E-mail: [e-mail] ^IIUniversidade de Brasília, Departamento de Engenharia Mecânica, 70910900 Brasília, DF. Brazil

^IIIYale University, Department of Chemical Engineering. USA

ABSTRACT

In this paper, three-dimensional boundary integral computer simulations of emulsions in shear flows are described. Results for ordered BCC emulsions with dispersed-phase volume fractions below the critical concentration are presented. Complex rheological features including: shear-thinning viscosities, normal stress differences, and a nonlinear frequency response are also explored. For deformable drops, pairwise collision produces net cross-flow displacements that govern self-diffusion of drops. We compute trajectories of two interacting drops in shearing and present interesting numerical simulations of three dimensional gravity-induced motion of two drops. The results also demonstrate the feasibility of simulating high-volume-fraction emulsions and foam.

Keywords: Boundary integral, emulsion, simulation, rheology, shear flow

Introduction

An emulsion is an immiscible mixture of two fluids, one of which is dispersed in the continuous phase, typically made by rupturing droplets down to colloidal sizes through mixing. The dynamics of drop deformation when an emulsion flows in low particle Reynolds number is of interest in a wide variety of fields including chemical and petroleum engineering. Furthermore, the behavior of drops in emulsions can be regarded as a prototypical motion of complex particles and capsules such as red blood cells.

There are several disadvantages which arise in connection with the application of existing non-linear continuum theories to actual flows problems (e.g. Noll, 1958 and Truesdell and Toupin, 1960). In the first place, the resulting constitutive equations contain numerous coefficients which, in general, cannot be measured with current experimental techniques. In addition, these coefficients, having been generated from a formal phenomenological approach, are usually devoid of any physical significance, and hence their dependence on the physical properties of the material is unknown. And finally, a dilemma is at once encountered when an attempt is made to accurately describe observations by retaining as much generality as possible in the existing theories, in that increased generality usually renders the problem unsolvable whereas too much simplification inevitably results in a poor description of the observed phenomena. It is clear that the purely mathematical approach to the formulation of constitutive equations should be supplemented by techniques which would result in simple, appropriate models for particular classes materials. In this paper this possibility is explored for one such class, consisting of an emulsion of droplets suspended in an incompressible Newtonian fluid. The rheology of this system is studied theoretically by a boundary integral method.

The boundary integral method relates velocities at a point within the fluid to the velocity and stress on the bounding surfaces. It is an ideal method for characterizing emulsion rheology since this problem is a direct consequence of drop deformation and drop orientation by the flow. This technique includes the reduction of problem dimensionality, the direct calculation of interfacial velocity, the ability to track large surface drop deformations, and the potential for easily incorporating interfacial tension (see Rallison, 1984 and Stone, 1994). The boundary integral formulation for creeping flows was theoretically described by Ladyzheskaya (1969) within the framework of hydrodynamic potentials. This integral equation method was first implemented numerically by Youngren and Acrivos (1976) in the study of Stokes flows past a particle of arbitrary shape. In addition, Rallison and Acrivos (1978) and Rallison (1981) performed numerical simulations to investigate deformation of viscous drops and bubbles in shear flows. Many aspects related to integral equation methods for free-boundary problems are described in a book by Pozrikidis (1992).

The principal difficulty with solving drop motion is that the position of drop interface is unknown a priori and must be determined as part of solution. Thus, the problem of determining the time-dependent drop shape is inherently non-linear. Emulsion rheology is difficult to predict or control due to the complex interplay between the detailed drop-level microphysical evolution and the macroscopic flow. There has been little progress toward a fundamental basis for understanding and predicting the rheology of such complex systems. Computer simulations provide a potentially valuable tool for helping to understand the microstructural mechanisms of emulsions flows. Although most work has been concerned with axisymmetric or two-dimensional interface drop deformations, which require numerical treatment of line integrals, several studies have considered the more difficult case of three-dimensional drop distortion (e.g. Mo and Sangani, 1994; Li et al. 1995; Loewenberg and Hinch, 1996, 1997; Cristini et al. 1998, Zinchenko et al. 1997; Cunha et al. 1999b). In particular, the work of Loewenberg and Hinch (1996) has provided a first basis for describing emulsion rheology by numerical simulations. The shear modulus, large-deformation elastic behavior, and quasistatic response in simple shearing flow have been studied for concentrated emulsions and foams by Kraynik and coworkers (e.g. Kraynik, 1988, Kraynik, Reinelt and Princen, 1991). Numerical simulations of concentrated emulsions are needed to help the interpretation of experimental observations (Mason et al. 1997) and to predict the rheology and microstructure of such complex flows.

The aim of this article is to demonstrate the application of boundary integral simulations in order to predict the microstructure effects on the rheology of emulsion flows. Using our simulation the linear viscoelastic response is computed for 3D ordered BCC (body centered cubic) concentrated emulsions for drop concentration f below the critical concentration f_c (maximum packing of spherical drops). For deformable drops, pairwise collision produces net flow cross-flow displacements that produces dispersion of drops in an emulsion. We compute trajectories of interacting drops in shearing emulsions and the motion of two interacting drops under gravity action. A substantial progress involving uniform expansion of a concentrated emulsion on a BCC lattice is also presented. Hydrodynamics interactions between neighbors drops cause drop distortion and large volumetric expansion rate of the drops promote the formation of thick liquid films.

Problem Formulation

Consider the creeping flow motion of an emulsion composed of suspended droplets with underformed radius a, undergoing oscillatory shearing flow field g_ocos(wt). Let g_o be the magnitude of the oscillatory shear (i.e. the shear rate), µ being the continuous-phase viscosity, and s the interfacial tension as shown in Fig. 1. All velocities are non-dimensionalized by s/µ and all lengths by a. The relevant parameters of the problem include the dispersed-phase volume fraction f, the dispersed- to continuous-phase viscosity ratio l, the capillary number Ca, the Bond number Bo and the Strouhal number W (i.e. dimensionless frequency). The capillary number is defined as the shear-rate normalized by the drop relaxation rate s/µ a

Bond number represents typical hydrostatic pressure variations relative to interfacial tension stresses: Bo = Dr g a² / s , where Dr is the density difference between drop and surrounding fluid. In the absence of buoyancy-driven motions of the drop Bo=0 (i.e. neutral buoyancy drops). The Strouhal number is relevant in the problems involving oscillatory shear flow. It is defined as being the ratio of the imposed oscillation frequency w to the drop relaxation rate:

Microscale Motion

The microhydrodynamics of the flow corresponding to drop length scale is treated at low particle Reynolds number, i.e.

Re = a²g_o / n « 1, where n is the kinematic viscosity of the fluid phase. Incompressible fluid motions at low Reynolds number are governed by the Stokes and continuity equations, namely

Here, T is the local Newtonian stress tensor given by

where the symbol Ñ=(¶/¶x_k)e_k denotes the partial differential operator nabla for a specified value of k (=1,2,3), t indicates transpose of the tensor Ñu, u is the Eulerian velocity field, I is the identity tensor, P=p-rg· x and P'=p-(r+Dr)g· x are the modified pressures of the continuous and dispersed phases respectively, and µ, r and l µ, r+Dr are the respective viscosities and densities of the phases. Dr is the density contrast between the two phases, g is the gravitational acceleration, x is the position vector and p denotes the mechanical pressure. Although time-dependence does not appear explicitly in Stokes equations, it is consistent to study time-dependent drop deformation. This quasi-steady assumption requires that a²/n « (1+l) µ a / s , where a²/n is a time scale for vorticity diffusion, and (1+l) µ a/s is the drop relaxation relaxation time. Physically, the quasi-steady approximation means that the fluid immediately adjusts to changes in the boundary location owing to rapid vorticity diffusion.

Boundary Conditions

The boundary conditions at a drop interface S with surface tension s require a continuous velocity and a balance between the net surface traction and surface tension forces (Pozrikidis, 1992). Mathematically, these conditions are expressed as following

where u^¥ is the imposed shear flow on the emulsion.

The traction jump Dt =[ n· T] at the interface is written as (Pozrikidis, 1992)

Here (I-nn) · Ñ denotes the gradient operator Ñ^s tangente to the inerface and n is the unit normal vector to S. Ñ · n defines the interface mean curvature 2k , where k is expressed as the sum of the inverse principal radii of curvature ½ (R₁^{-1 +} R₂^-1 ).

The force balance at the interface, Eq.(7), includes both a normal traction jump (n·Dt) n and a tangential stress jump (I-nn) ·Dt due to gradients of surface tension along the interface. The dynamical effect due to interfacial tension gradients are known as Marangoni motions, caused by the presence of surfactants in the fluid and/or temperature variations . In real flow systems, one must often expect surface tension to vary from point to point on the interface. Consequently, it is important to consider how gradients of s may influence drop deformation. To describe Marangoni stress effects, an additional convection-diffusion equation is necessary in order to determine the surfactant distribution along the interface (Stone, 1994). The formulation herein, however, considers the case of clean interface deformable drops. Under this condition the interface is treated as having uniform surface tension. The traction disconinuity given in Eq.(7) reduces then to a contribution from capilary pressure and gravity only:

An additional kinematic constraint is necessary to relate changes in the interface position x_i to the local velocity u (x_i). Therefore the interface evolution is described in a Lagrangian representation as

where M is the number of interfacial marker points on the drop surface.

Macroscopic Motion

It is assumed that the macroscopic length is much larger than the drop size a. Thus, the emulsion may be seen as a homogeneous continuous fluid. Accordingly, the macroscopic flow U field is governed by the Cauchy and continuity equations respectively

Here, is the effective density of the emulsion, S is the emulsion volume-averaged stress tensor which is computed by direct numerical simulation of the detailed microstructure. Although the ambient fluid and drop phases are assumed Newtonian, S is generally non-Newtonian as a result of drop deformation and drop orientation in the flow direction. The effective stress tensor S of an emulsion is given by the average of the actual stress T over a volume V of the emulsion (Landau and Lifshitz, 1987):

Following Batchelor (1970) and Hinch (1977), if the suspending fluid is Newtonian and the flow occurs at low particle Reynolds numbers, then

The first two terms on the right-hand side are the contributions from the background continuous phase (i.e. a Newtonian effect) with an average pressure á Pñ and an average rate of strain áDñ. The average shear rate of the flow is defined as . Here, tr denotes the trace of a tensor, n is the average number of drops per unit volume (N/V) and E is the average stress contribution of the dispersed phase, and it value depends on the size, shape, orientation the drop, and on the positions, sizes, shapes and orientations of the other dops in the emulsion (i.e. a configuration of the surrounding drops):

The subscript m denotes a particular drop of surface S_m , and the summation in Eq. (13) is calculate for m varying from 1 to N.

In simple shear flows, namely, u^¥ = (g_ox₂,0,0), the continuous phase contributes to the emulsion average shear stress as

An effective shear viscosity µ^* = S₁₂ / g_o is defined as

The last term on the right-hand side of Eq. (15) corresponds to the shear thinning effect in the emulsion due to drop deformation. Nonzero first and second normal stress differences result entirely from the drop contribution

Linear Viscoelastic Response of Concentrated Emulsions

Emulsions posses microscopic mechanisms for both elastic and viscous dissipation. The energy storage and dissipation per unit of volume can be represented by a frequency-dependent complex viscosity h^* (f,w) =h' (f,w) - i h'' (f,w), which is defined only for strain amplitude sufficiently small so that the shear stress is linear in the rate of strain. Here f is the droplet volume fraction and w is the frequency. The real part h' (f,w) = w G^ll (f,w) is the in-phase ratio of stress with respect to an oscillatory rate of strain, and reflects dissipative mechanisms, whereas the imaginary part

h' (f,w)=w G^l (f,w), is the out-of-phase ratio of stress and reflects elastic mechanisms. G^l (f,w) and G^ll (f,w) are called storage modulus and loss modulus respectively. In this work we explore the behavior of these quantities over a wide range of f and w in order to characterize the importance of the dissipative and elastic mechanism as the droplets become packed and deformed.

The average flow shear-rate for the case of an oscillatory shear is given by

where t denotes time, and w and g_o the frequency and the amplitude of the oscillations, respectively. Then, for an emulsion undergoing a simple oscillatory shear the far field flow will be given by u¥ = (g_o x₂ exp(iwt),0,0).

When the rate of strain g_o in the system is small, the stress response is linear. The stress oscillates with the same frequency but not necessarily in phase with the forcing. Based on a linear viscoelastic formulation for dilute limit of an emulsion Cunha et al. (1999a,b) developed a model of a single relaxation time. According to this model, the in-phase part of the response (viscous dissipation) is found to be

whereas the out-of-phase ratio of stress reflecting elastic effect is given by

where h_o =(1+5l/2) / (1+l) is the viscosity in the limit of low frequency (Taylor viscosity, Taylor, 1932), whereas h_¥= 5(l-1)/ (2l+3) is the viscosity in the limit of high frequency, corresponding to the viscosity of an emulsion composed of spherical blobs (surface tension is not an important effect in this case). In Eq. (18) and Eq. (19), t_o = 40(2l+3) (19l+16) / (l+1) is a single drop relaxation time.

Following the above model, the drop stress contribution, for higher concentrations, can be expressed in an appropriate form as a purely viscous contribution plus a complex stress contribution from the surface tension, G.

A complex viscosity h^* is defined from G = h^* (f,w) g (t), and G is obtained from the stress relation function spectrum F(f,t) by using the convolution integral

Alternatively one may define the complex viscosity in terms of the stress relaxation function, directly. By Fourier transformation it possible to show that

Now, the average relaxation time is given by the following Integral

It is instructive to note that, for dilute limit of emulsions, =t_o and F(0) = F_o= h_o/t_o. The simplest viscoelastic model studied is the Maxwell model of one relaxation time. Namely,

Dilute emulsion are well described by such a model. On the other hand, it will be shown in the result section that a multiple Maxwell relaxation time model (t_k, F_k ; k=1,2,...), as defined below, is needed in order to describe the rheological behavior of concentrated emulsions undergoing oscillatory shear with small amplitudes and arbitrary frequencies.

Description of the Simulations

All simulations rely on the boundary integral method. The initial interfaces of the drops are triangulated by subdividing the edges of an icosahedron into f equal segments, thereby forming subtriangles on each face of an icosahedron as shown in Fig. 2 (a,^b). The additional vertices are displaced radially outward onto the spherical surface that circumscribes the icosahedron. By this procedure, the drop surface is discretized into nt=20 f² flat triangles of equal area. All of the 10f²+2 vertices have coordination number 6 except for the evenly spaced icosahedron vertices that have coordination number 5 (Almeida and Cunha, 1998).

We have developed three-dimensional boundary integral simulations that are capable to describing the dynamics of concentrated emulsions in oscillatory shear. Accordingly, the evolution of N deformable drops is described by time-integrating the fluid velocity u(x_o) on a set of interfacial marker points x_o on each drop surface. The fluid velocity is governed by the second-kind integral equation on the interfaces S_m (m=1,...,N) of all drops (Rallison and Acrivos, 1978)

with

where b= 1/(l+1), y_o is the coordinate in the direction e₂ of the velocity gradient for the prescribed shear flow, x Î S_m is the integration point, and r = x-x_o . Equation (26) has been made dimensionless using the underformed drop size a for characteristic length and the cappilary relaxation velocity s/µ. Ca^*(t)=Ca cosWt for oscillatory shear, and Ca^*=Ca for stationary shear. The dimensionless surface traction is

For isolated drops in Stokes flow, the kernel functions J e G in Eqs. (26) and (27) are the free-space Green's functions called stokeslet and stresslet kernel function in the context of creeping flow hydrodynamic:

For concentrated emulsions, however, periodic boundary conditions are imposed summing the free-space kernel functions on a lattice using Ewald-summation given in Beenaker (1986). The resulting periodic Greens' functions G^P and J^P are discussed in details by Loewenberg and Hinch (1996) and in a recent work of Cunha, Abade, Sousa e Hinch (2002) for the threedimesional case. These periodic functions are obtained by Ewald summation using accurate computationally-efficient tabulation of the nonsingular background contribution. A strained cubic lattice is used to describe shear flow; by symmetry only strains in the range [-1/2,1/2] are needed.

The boundary integral simulation here uses an adaptive drop interface discretization which depends on the instantaneous drop shapes only. It is independent of fluid velocity and also of the drop deformation history. In particular, the algorithm incorporates a prescribed marker-point density function (interface area per marker point) that is used to define an equilibrium edge length between marker points. Following Cristini et al. (1998), we define a local desired density on the drop surface as being r_M = Ck². The constant C is related to resolution of a spherical shape. Let N₀ be the number of marker points chosen to discretize a spherical drop of radius a, then C~ N₀/4p. This density function resolves the minimum local length scale everywhere on the drop interface. Equilibration velocities for marker points are determined by the resultant of local spring-like tensions projected onto the drop interface. The resulting dynamical system of damped massless springs has a well-defined minimum energy equilibrium state that is attained in several iterations after every fluid dynamic displacement of the marker points. Marker point equilibration is an inexpensive O(N) calculation compared to the time-controlling O(N²) fluid velocity calculation. The interface discretization algorithm maintains optimal marker-point connectivity using a local reconnection rule. Accordingly, the edge between adjacent triangles is switched to connect opposite vertices, if the resulting edge is shorter. By maintaining a nearly equilateral triangulation, this procedure maximizes the time step for the time-controlling fluid dynamical calculations. The normal vector and curvature were calculated by the local surface-fitting algorithm of Zinchenko et al. (1997).

The fluid velocity on the drop interfaces is obtained by an iterative solution of Eq. (26) using the GMRES algorithm (a generalization of the conjugate gradient method for non-symmetric matrices) to achieve convergence for the closely-spaced interface configurations that characterize dense emulsions. Evolution of the deformable drop microstructure and associated stresses is obtained by time-integrating the fluid velocity on the interfacial marker points dx_k/dt, k=1,N. The positions of marker points are evolved by a second order Runge-Kutta scheme. Following Rallison (1981) an appropriate time step that is proportional to the shortest edge length is set to ensure stability.

Equation (26) has eigensolutions that cause unphysical changes in the dispersed-phase volume at small viscosity ratios, corrupt numerical solutions at large viscosity ratios, and slow the iterative convergence. These effects were eliminated by implementing Wielandt eigenvalue deflation described by Pozrikidis (1992) in order to purge the solutions corresponding to l=0 and l ® ¥. The resulting surface integration algorithm is economical and (1/N) accurate, consistent with the triangulated representation of the drop interface. For efficient trapezoidal-rule integration the Eq. (26) was transformed into a singularity-subtracted form described by Loewenberg and Hinch (1996) and solved by Jacobi iteration.

Finally, the adaptive surface triangulation algorithm, described by Cristini et al. (1998), was extended to construct efficient simulations of dense emulsions or foam. Accordingly, a new marker point density function was defined that resolves the minimum local length scale everywhere on the drop interfaces. For the proposed problems here the minimum local length scale may depend on the local curvature or local film thickness h. Thus, the marker point density function should be generalized to

where R₁, R₂ are the local principal radii of curvature, and C₁ e C₂ are O(1) constants whose precise value is unimportant. The proposed marker point density function resolves the rim of dimpled regions where the film thickness varies rapidly and the lubrication length scale (h R)^1/2 in regions where the film thickness is slowly-varying (R=min[R₁,R₂]). Only rim regions, not flat regions, require high resolution. An inspection of the result in Fig. 3 which was obtained using the new density function, illustrates this point: high resolution is needed on the plateau borders and junctions, not on the flat films between drops. This observation is important for ensuring the feasibility of accurately resolving the dynamics of the closely-spaced interfaces that characterize the dense systems that we have explored when simulating emulsion expansion. The grid remains well-structured during entire simulations. Nearly-equilateral triangles are maintained by edge reconnection and local equilibration. As the drop deformation increases marker points are continuously added or subtracted to guarantee constant interface resolution according to the density function defined in Eq. (30). For simulations of emulsion expansions Cunha, Sousa e Loewenberg (2003) have recently developed a new boundary integral for compressible Stokes flow finding the appropriate function F(x_o) for this case as being

where the viscous free-surface flow problem of emulsion densification should include the new parameter L= µ D a / s defined as the capillary expansion of the drops, and D is the volumetric expansion rate of the emulsion. Here, f(f )=1+l(f ^-1-1).

Results and Discussions

In this section we show simulation results of the deformation of two interacting drops in shear flow, the deformation of two buoyant drops, emulsion rheology and emulsion expansion for large drop concentration.

Drop Deformation

In Fig. 4 stationary drop shapes are depicted for steady shear and different values of drop viscosity, l=3.8, 5.0, 10. An inspection of the drops illustrates that the adaptive discretization resolves regions of high curvatures where a high density of marker points is needed. In particular, Fig. 4 shows that drops with sufficiently high viscosity, the characteristic time for drop deformation, (1+l) µ a/s, exceeds the characteristic time for drop rotation, 1/g _o; hence the drop rotates from the extensional to the compressional quadrant of the velocity gradient before being significantly deformed, and this prevents breakup. Cunha, Toledo and Loewenberg (1999a) have shown theoretically that dilute emulsions with high viscosity drops are shear thinning and have normal stresses at high shear rate. This behavior being a direct consequence of drop deformation and alignment with the flow shown in Fig 4. It is seen that an elongation of drops in the direction of the flow offers less resistance for shearing the emulsion and produce normal stress differences.

Rheology of a High Viscosity Emulsion in Steady Shear

The volume-averaged stress tensor in shear flow can be expressed as the stress that would result from only continuous-phase fluid and the extra stress E that results from the dispersed-phase drops as given in Eqs. (12) and (13). Locally small drop Reynolds numbers, neutrally buoyancy (Bo=0), and negligible Marangoni stresses are assumed. Under these conditions the traction jump (normalized by s / a) in Eq. (28) reduces simply to Dt = 2k (x) n. The surface integral (13) requires only the shape and velocity distribution on the drop interfaces which is directly obtained by solution of Eq. (26) from the boundary integral formulation described above. Therefore, the main part of the bulk stress of an emulsion in simple shearing motion is the calculation of the stress contribution of the drops E due to the dipole stresslet that each drop torque free generates in the ambient fluid as it flows.

Figure 5 shows a comparison between numerical results of the drop stress contribution from the boundary integral simulation and a high l drop theory developed in our previous work (Cunha, Toledo and Loewenberg, 1999a). The shear thinning and the elastic behaviors (i.e. normal stress differences) predicted by the high viscosity drop theory for l=10 and for l=20 agree well with the boundary integral simulations. For l=5, however, it is observed a significant discrepancy between the results, because the drop is not slightly deformed as assumed in the theory. Fig. 5a shows asymptotic limits at small and high dimensionless shear rate (i.e. Ca), and a shear thinning region where the effective viscosity S₁₂ / g depends on the shear rate. At low shear rate the emulsion behaves like a Newtonian fluid with Taylor viscosity (1932) (i.e. spherical drops with high surface tension). In the limit of high shear rate the emulsion is also Newtonian but with an effective viscosity due to blob particles . This limit goes to zero due to the subtraction of the blob viscosity from the drop stress contribution. In the case of spherical blobs the surface tension does not exert any influence on the emulsion nonlinear response.

The last two plots in Figs. 5b and 5c show the elastic behavior of the emulsion in terms of the first and second normal stress differences N₁ and N_2, respectively. It is seen that for all shear rates, N₁ while is a positive function of the Ca, N₂ is less than (e.g. N₂@ 1/4 N₁, for l=5) and always negative; a typical behavior of elastic liquids (Bird et al. 1987). It can be seen that at low shear rates the normal stress differences appear as a second order effect (nonlinear), proportional to Ca², while S12 scales like Ca. The simulations demonstrates that drop deformation and the rapid alignment with the flow for high viscosity drop emulsion (as illustrated in Fig. 4) lowers the shear stress and produces elastic stresses. The reduced collison cross-section imposed by drop deformation allows them to slide past each other with less resistance and the orientation of the drops in the flow direction produces stress anisotropy and a large normal stress differences. These qualitative features are consistent with stress contributions for high viscosity drops theory (l ® ¥) in stationary shear given by the asymptotic limits of small capillary number and dilute emulsion where: S12=fCa (5/2-3/2l), N₁=10Ca²/lc² and N₂=- (29/14l)Ca²/c with c=20/19 (Schowalter, Chaffey and Brenner 1968, Cunha et al. 1999a).

Rheology of a High Viscosity Emulsion in Oscillatory Shear

Figure 6a shows a typical numerical result of an isolated high viscosity drop (l=5) deforming in oscillatory shear with a dimensionless shear rate (capillary number) Ca=3 and W=0.5. The motion is periodic and the drop deforms, rotates, but it does not break under the conditions of this simulation. On the other hand the result depicted in Fig. 6 (b) indicates that drops with smaller viscosity, l=1, can break in oscillatory shear since the strain amplitudes are large. We have explored the limiting high shear-rate rheology of diluted emulsions with high viscosity drops. Even at high capillary number (i.e. strong flow) no large deformation results for high viscosity drop, because the extended part of the shape is spun round into the compression before it has extended far. When the equilibrium is such that the slightly deformed droplet has its principal extended axis in the direction of the flow, the viscous dissipation within the external fluid is decreased as a direct consequence of less distortion in the flow streamlines, and a shear thinning behavior may be observed with the presence of the first and the second differences of normal stresses.

Figure 7 shows a complex rheology corresponding to a nonlinear frequency response for high strain amplitudes (Ca/W=50). The horizontal dot lines represent the steady value of I₁ and I₂. It was found that the behavior of high viscosity drops depends on three rates: the shear-rate g, the oscillation frequency w, and the characteristic time for drop deformation that is (1+l)µa/s. At low shear-rates, the system exhibits a linear viscoelastic response and a nonlinear behavior occurs for moderate and large shear-rates. In this plot the analytical formulae of the drop shear stress contribution and the first normal stress difference for a dilute emulsion of high viscosity drops (l ® ¥) in oscillatory shear were developed by Cunha et al. (1999a, b):

where

The integrals I₁ and I₂ were evaluated numerically using a trapezoidal rule. The spacing of 0.01 between points produced an error less than 0.1%. In the asymptotic limit of steady shear (W ® 0) we find the following expressions:

Simulations of Two Drops Interacting in Simple Shear

Our numerical method allows detailed computation of interactions between deformable drops. In this part of the simulations we show simulation results for typical collision of drops in shearing motions. Actually, the motion of interacting drops in shear flows may be regarded as a first prototypical motion of red cells in microvessels in which Reynolds is about 0.01 and two cells interactions dominate.

In the absence of Brownian motion, inertia and interparticle forces, two equal smooth spheres collide in simple shear flow in a reversible way returning to their initial streamlines. On the other hand experiments carried out by Arp and Mason (1977) , theory (Cunha and Hinch, 1996) and numerical simulations (Loewenberg and Hinch 1996) show that in suspensions of rough rigid particles or in emulsions of deformable drops, the net particle displacements after interacting may produce a source of mixing which may be characterized by a self-dispersion process. In order to check this phenomenon we have computed pair interaction between two equal deformable drops undergoing simple shear flow for Ca=0.2, l=1.0, and initial horizontal and vertical offsets x₀=-10 and y=0.4. A sequence of six stages of the colliding drops at successive times is illustrated in Fig. 8. The net cross-flow displacement resulting from the drop interaction is apparent. In fact, it demonstrates that interactions between two deformable cells tend to increase the cross-flow separation of their center.

In Fig. 9 the relative trajectory is plotted to capture the approach, collision and separation between the drops as shown in Fig. 8. It is observed that the drops separate on streamlines further apart on their approach with a net transversal displacement Dy @ 0.9. Thus, the interaction is irreversible, and repeated collisions led to increasing values of the net displacement until drop interactions became negligible. This prediction may have direct relevance to the accumulation of platelets in a vessel wall due to shear induced migration. These result are in good agreement with the experimental observation of Guido and Simeone (1998), who investigated binary collision of drops by computer-assisted video optical microscopy and image analysis techniques.

Results are also obtained for the deformation parameter D, defined as (L-B) / (L+B), where L and B are, respectively the major and minor axis of the deformed drop during collision. According to Fig. 10, D is constant during the approach (configuration a) and, the deformation of the drops are maximal when the drops are pressed together along the compressional axis of shear (configuration b) and minimal when they are drawn apart along the extensional axis (configuration c). The extensional action of the surrounding fluid leads to a second increase of deformation up to the maximum when the drops separate (configuration d). The two drops then relax to the steady shape in the absence of interaction. Eventually D reaches a steady state at the same plateau value as before the collision. The shear stress contribution of the drops follows the similar behavior of deformation, whereas elastic stress contributions are maximal when the drops are aligned with the flow.

Relative Gravity-Induced Motion of Two Drops

The simulation presented in Fig. 11 shows the relative gravity-induced motion of two drops at small Reynolds number for l=1 and Bond number, Bo=1. The particle initial relative position was (1.2,1,0) and the radius of the large drop 1.5a. The results indicate a significant drop deformation during the interactions, demonstrating the possibility of a low viscosity drop breakup as a result of pair interaction in buoyancy-driven motion. As shown in the frames of Fig. 12, lubrication stresses between the deformed drops prevent coalescence; the smaller drop slides when it is passing the larger one. The last frame shows the smaller drop stretched in the straining flow behind the larger one; the small drop is swept around the large drop. Physically, the flow created by the leading drop (larger drop) tends to flatten the trailing drop into an oblate shape. A simple scaling arguments gives at leading order the total strain experimented by the trailing drop as being

Bo a₁ / d, where a₁ is the radius of the leading drop and d is the separation distance. Therefore for separation distances of 5 drop radii and Bo=1, we can expect strains of 0.2 a₁. Drop coalescence is not observed in the simulations since modelling of the interfacial phenomena, such as van der Waals attraction, is not included in the present studies. These results are in qualitative agreement with the experimental sequence of centimeter-size bubbles interacting in corn syrup carried out by Manga and Stone (1993) and with numerical simulations Zinchenko et al. (1997) , who has also investigated motion of drops in sedimenting emulsions.

Rheology of Concentrated Emulsions

Next, results of 3D step strain simulations for drop concentration below the maximum packing of spherical drops (f < f_c=p 3^1/2 / 8 » 0.68) and viscosity ratio l=1 are presented. The emulsion was subject to a step strain g(t)= g₀d (t), allowed to relax and the response was measured. Figure 12 shows the dimensionless stress relaxation function as a function of the time normalized by the average relaxation time for three different concentrations: dilute (f<0.3), 0.5 and 0.66. It is seen that the dilute system has a simple exponential decay and also that one Maxwell relaxation time model describes well the response of the system. For all simulated volume fraction it was observed a long time exponential decay. The short time scales, on the other hand, show some multiple time scales, especially at the highest simulated volume fractions (i.e. f=0.66).

A typical behavior of the average relaxation time is shown in Fig. 13. The average relaxation time is normalized by the drop relaxation time and it is plotted as a function of the dimensionless volume fraction . It is seen that the average relaxation time diverges as f is going close to the critical concentration. This divergence is closely related to the divergence of the low frequency viscosity when f ® f_c, and it is described by the lubrication resistance between the drops. In the regime of high concentrations, the interactions between the drops are dominated by high order multipole and lubrication forces acting between them. Consequently, a Maxwell model of one relaxation time, in this case, should not be applied for describing the emulsion behavior.

Relaxation times can be extracted directly from our numerical simulations. The spectrum of these relaxation times is shown in Fig. 14. It is seen that the long time mode diverges and dominates the average relaxation time. At higher concentrations it is observed the presence of multiple modes (i.e. multiple time scales) which may be attributed to the higher order multipole interactions between the drops and the lubrication resistance between the drops. In addition , Fig. 14 indicates that the frequency response of a dilute emulsion is very well-described by one relaxation time, whereas four relaxation times are needed to describe the system response at 66%. The results obtained for the real (i.e. viscous contribution) and imaginary (i.e. elastic contribution) parts of the complex viscosity are presented in Fig.15. It is seen that the multiple relaxation time observed in the spectra shown in Fig. 14 seems to be sufficient for describing the viscoelastic response of a BCC emulsion at 66% drop volume fraction. Thus, a Maxweel model of four relaxation times fits very well the numerical results.

At all concentration simulated , we estimate the low-frequency regime in which G^l (f,w) depends slightly on the frequency. The low frequency drop presumably reflects the very slow relaxation of the elastic structure of the emulsion, while the high-frequency limit reflects the fact that the emulsion is comprised solely of fluids, whose viscous behavior dominates and surface tension is unimportant (i.e. the drops behave as blobs). In order to predict an equivalent static shear modulus of the G_o^l (f) emulsion we use the low-frequency limit given by the simulations. From the results shown in Fig. 15, G_o^l (0.66) » 2 x 10^-3s/a (units of s/a). The simulations capture a rich rheological behavior of emulsion flows; a yield stress can be recovered in the zero-frequency limit that is in qualitative agreement with the experimental value of Masson et al. (1997) at the same concentration.

Emulsion Expansion

Emulsion densification with prescribed nonequilibrium osmotic pressure has also been explored since it corresponds most closely to the experimental procedure that has been developed by Mason et al. (1997). Under these conditions, the drop expansion proceeds until the applied osmotic pressure is balanced by the equilibrium osmotic pressure of the emulsion. We have made substantial progress on this problem by using the boundary integral method described earlier. Thus far, we have only considered the case l=1 because the numerical implementation is simpler. Results of the present simulations for L =0.5 are contained in Fig. 16a,^b, which illustrates the final structure generated by the expansion of drop shape from spheres to Kelvin cells and a Weaire-Phelan microstructure, respectively. The spherical drops on a BBC lattice make first contact with eight nearest neighbors and form precursors of hexagonal faces when f>f_c=0.68 (i.e. maximum packing for a BCC structure). Once the emulsion expansion has stopped, the drop shape continues to relax toward equilibrium; the emulsion expand and then relaxes for twenty time units promoting the formation of thick liquid films and a redistribution of liquid from the gap to the plateau border relaxation.

Conclusions

In this work we have presented numerical simulation results for the rheology of emulsions under steady and unsteady shear flows. The results revealed a complex rheology of an emulsion, including shear thinning and normal stresses, at dilute and moderate regime of drop volume fraction below the maximum packing configuration. The main conclusions of this work are: i) The results have demonstrated that drop deformation and alignment with the flow lowers the shear stress and produce large normal stresses. ii) For an emulsion of high viscosity drops the characteristic time for drop deformation exceeds the typical time for drop rotation; hence, the drop rotates from the extensional to compressional quadrant of the shear before being significantly deformed. iii) By tracking the relative motion of two colliding drops in a dilute emulsion, it was found that the distance between the centers of mass along the velocity gradient did not restore, after separation, to the value before collision, but became larger. It suggested that drop deformation during collision provide the mechanism for diffusion of drops across streamlines firstly observed by Cunha and Hinch (1996) for rough particles and Loewenberg and Hinch (1997) for deformable drops. iv) The buoyancy-driven interactions between two drops have shown that for Bo=1 and a moderate viscosity ratio l, the small drop just swept around the larger drop without coalescence. v) The stress relaxation function decays exponentially at long time for all tested cases; iv) The time relaxation spectra were found to be discrete for ordered BBC emulsion; vi) The average relaxation time diverges and this divergence is related to the lubrication resistance between the drops when f ® f_c; vi) Dilute emulsion f< 0.3 are well described with a Maxwell model of one relaxation time, whereas a superposition of four relaxation times is needed to describe successfully the viscoelastic response of an emulsion with drop volume fraction around 66%. The bulk elastic property of the emulsion was probably the most important feature of our results. In the present simulations , such elastic behavior resulted of course from the presence of the finite interfacial tension which always acts so as to oppose any deformation of a drop from the spherical shape. The emulsion was seen to exhibit properties usually associated with viscoelastic fluids. This result is certainly of interest since both of the component fluids are Newtonian liquids. vii) The numerical approach was also applied for generating a dense compressed emulsion structure in viscous flows with periodic boundary conditions. The results demonstrate the feasibility of simulating high-volume-fraction systems. Results from investigation on the rheology of disordered emulsions for concentrations above the critical volume fraction will be the subject of future publications.

Acknowledgments

The authors are grateful to the CAPES-Brazil, CNPq-Brazil and CT-PETRO-FINEP for their generous support of this work. We would like to thanks Dr. Vittorio Cristini and Dr. Jerzy Blawzdziewicz for helpful discussions on their algorithm for adaptive mesh restructuring of three-dimensional drop surface.

Paper accepted November, 2002. Technical Editor: Aristeu da Silveira Neto

• Almeida, M.H.P. and Cunha, F.R., 1998, "Numerical simulation of noncolloidal deformable particles at low Reynolds number", Proceedings of the Brazilian Congress on Engineering and Thermal Sciences, Vol. 11, pp. 963-970.
• Arp, P.A, and S.G. Mason, 1977, "The kinetics of flowing dispersion, Doublets of rigid spheres", J. Colloidal Interface Sci., Vol. 61, pp. 44-61
• Batchelor, G.K., 1970, "The stress in a suspension of force-free particles", J. Fluid Mech., Vol. 41, pp. 545-570.
• Beenaker, C.W.J., 1986, "Ewald sum of the Rotne-Prager tensor", J. Chem. Phys., Vol. 85, pp. 1581-1582.
• Bird, R.B., Armstrong, R.C., Hassager, O., 1987, "Dynamics of polymeric liquids" John Wiley & Sons, New York.
• Cristini, V., Blawzdziewicz, J. and Loewenberg, M., 1998, "Drop breakup in three-dimensional viscous flows", Phys. of Fluids, Vol. 10, pp. 1781-1783.
• Cunha, F.R. and Hinch, E.J., 1996, " Shear-induced dispersion in a dilute suspension of rough spheres, J. Fluid Mech., Vol. 309, pp. 211-223.
• Cunha, F.R., Toledo, R. and Loewenberg, M., 1999a "Hydrodynamic and Rheology of a Dilute Emulsion with High Viscosity Ratio", Proceedings of the Brazilian Congress of Mechanical Engineering", Vol. 1, pp. 1-10.
• Cunha, F.R., Neimer, M., Blawzdziewicz, J. and Loewenberg, M., 1999b, "Rheology of an emulsion in shear flow", Proceedings of the AIChE", Fundamental Research in Fluid Mechanics: Particulate and Multiphase Flow II, 124j, pp. 1-9.
• Cunha, F.R., Abade, G.C., Sousa, A.J., and Hinch., E.J, 2002, "Modeling and direct simulation of velocity fluctuations and particle-velocity correlations in sedimentation", ASME-J. Fluid Engn, Vol. 124, pp. 957-968..
• Cunha, F.R., Sousa, A.J., and Loewenberg, M., 2003, "A mathematical formulation of a boundary integral equation for compressible Stokes flows" Computational and Applied Mathematics (in press).
• Davis, R.H., Schonberg, J.A. and Rallison, J.M., 1989, The lubrication force between two viscous drops", Phys. Fluids, Vol. A 1, pp. 77-81.
• Guido, S. and M. Simeone, M., 1998, "Binary collision of drops in simple shear flow by computer-assisted video optical microscopy", J. Fluid Mech., Vol. 357, pp. 1-20.
• Hinch, E.J., 1977, "An averaged-equation approach to particle interactions in a fluid suspension, Vol. 83, pp. 695-720.
• Kraynik, A.M., 1988, "Foam flows", Ann. Rev. Fluid Mech., Vol. 20, pp. 325-357.
• Kraynik, A.M., Reinelt, D.A., and Procen, H.M., 1991, "The nonlinear elastic behavior of polydisperse hexagonal foams and concentrated emulsions," J. Rheology, Vol. 35, 1235-1253.
• Ladyzheskaya, O.A., 1969 "The mathematical theory of viscous incompressible flow", Gordon & Breach, New York.
• Landau, L.D and Lifshitz, E.M, 1987, "Fluid Mechanics, Course of Theoretical Physics", 2sd ed, Pergamon Press.
• Li, X., Zhou, H. and Pozrikidis, C., 1995, "A numerical study of the shearing motion of emulsions and foams", J. Fluid Mech. Vol. 286, pp. 379-404.
• Loewenberg, M. and Hinch, E.J., 1996 "Numerical simulations of a concentrated emulsion in shear flow", J. Fluid Mech., Vol. 321, pp. 395-419.
• Loewenberg, M. and Hinch, E.J., 1997, "Collision of deformable drops in shear-flow", J. Fluid Mech., Vol. 338, pp. 229-315.
• Manga, M. and Stone, H.A, 1993, "Buoyancy-driven interactions between two deformable viscous drops", J. Fluid Mech, Vol. 256, pp. 647-683.
• Mason, M.D., Lacasse, S.G. Grest, D. Levine, J. Bibette, D.A. Weitz, 1997, "Osmotic pressure and viscoelasticity shear moduli of concentrated emulsions", Physical Review E, Vol. 56} pp. 3150-3166.
• Mo, G. and A.S. Sangani, A.S., 1994, "Method for computing Stokes flow interactions among spherical objects and its application to suspension of drops and porous particles", Phys. Fluids, Vol. 6, pp. 1637-1652.
• Noll, W., 1958, "A mathematical theory of the mechanical behaviour of continuous media", Arch. Rat. Mech. Anal. Vol 2. pp. 198-226.
• Pozrikidis, C., 1992, "Boundary integral and singularity methods for linearized viscous flow", CUP, Cambridge.
• Rallison, J.M., and Acrivos A., 1978, "A numerical study of the deformation and burst of a viscous drop in an extensional flow", J. Fluid Mech., Vol. 89, pp. 191-200.
• Rallison, J.M., 1981, "A Numerical study of the deformation and burst of a viscous drop in general shear flows", J. Fluid Mech., Vol. 109, pp. 465-482.
• Rallison, J.M., 1984, "The deformation of small viscous drops and bubbles in shear flows", Ann. Rev. Fluid Mech., Vol. 16, pp.45-66.
• Schowalter, W.R., Chaffey, C.E. and Brenner, H., 1968, Rheological behavior of a dilute emulsion. J. Colloidal Int. Sci., Vol 26, pp. 152-160.
• Stone, H.A., 1994, "Dynamics of drop deformation and breakup in viscous fluids", Ann. Rev. Fluid Mech., Vol. 26, pp. 65-102.
• Taylor, G.I., 1932, "The viscosity of a fluid containing small drops of another fluid", Proc. R. Soc. A., Vol. 138 pp. 41-48.
• Truesdell, G.I., and Toupin, R.A., 1960, Handbuch der Physik, Vol 3., part 1, Springer, Berlin.
• Youngreen, G.K. and Acrivos, A., 1976, "On the shape of a gas bubble in a viscous extensional flow", J. Fluid Mech, Vol. 76, pp. 433-442
• Zinchenko, A.Z., Rother, M.A and Davis, R.H., 1997, "A novel boundary-integral algorithm for viscous interaction of deformable drops", Phys. Fluids, Vol. 9, pp. 1493-1511.

Publication Dates

• Publication in this collection
  18 Mar 2004
• Date of issue
  Mar 2003

History

• Received
  Nov 2002