A general approach for modelling complex fluids: its application to concentrated emulsions under shear

Colloids and Surfaces A: Physicochemical and Engineering Aspects 144 (1998) 305–318

A general approach for modelling complex fluids: its application to concentrated emulsions under shear Xue-Feng Yuan a,b,1*, Masao Doi b a Kao Institute for Fundamental Research, Kao Corporation, 1-3, Bunka 2-chome, Sumida-ku, Tokyo 131, Japan b Department of Applied Physics, Faculty of Engineering, University of Nagoya, Furo-cho, Chikusa, Nagoya 464-01, Japan Received 9 April 1998; accepted 1 July 1998

Abstract We present a general Lagrangian–Eulerian approach for modelling complex fluid flows. Our method can track the microstructure evolution of complex fluids under flow in a natural physical way. It is especially useful for mesoscopic modelling of complex fluids in which the characteristic time and length scales are usually several orders of magnitude higher than those for simple fluids. The method is illustrated by studying hydrodynamically interacting two-dimensional concentrated emulsions under Couette flow. We have carried out the simulations of concentrated emulsions under shear over a range of capillary number (0.1≤Ca≤0.6) with a fixed area fraction w=0.4 and viscosity ratio l=1. The results show that concentrated emulsions exhibit the characteristics of viscoelastic fluids: strong shear thinning of shear stress and non-zero of first normal stress difference. Deformable drops can glide past each other with less resistance than in the case of hard suspensions and hence reduce shear stress. We have found that the dependence of the first normal stress difference on shear rate changes from quadratic to linear as the shear rate increases. In the regime of larger capillary numbers, highly inhomogeneous behaviour of drops and the dumbbell shape of drops can be clearly observed. © 1998 Elsevier Science B.V. All rights reserved. Keywords: Emulsions; Couette flow; Voronoi; Rheology; Capillary number

1. Introduction Complex fluids, such as emulsions, suspensions, polymer mixtures and lyotropic liquid crystals, are multi-component heterogeneous systems. Modelling of complex fluid flows presents a considerable challenge. The characteristic time of a typical complex fluid is of the order of 1 s and its length scale is in a range of microns to millimetres. Restricted by current computer power, the applica* Corresponding author. E-mail: [email protected] 1 Present address: H.H. Wills Physics Laboratory, Royal Fort, Tyndall Avenue, Bristol BS8 1TL, UK.

tion of conventional Monte Carlo and molecular dynamics techniques for atomic or molecular fluids to complex fluids is still quite limited. On the other hand, finite-element and finite-difference methods are often adopted to solve continuum partial differential equations for engineering problems provided an appropriate constitutive equation for given materials and boundary conditions are known. Unfortunately, the assumption of a homogeneous fluid for continuum methods breaks down for many complex fluids in which their heterogeneous nature cannot be ignored over the time and length scales of their typical physical process. Furthermore, no constitutive equation is available

0927-7757/98/$ – see front matter © 1998 Elsevier Science B.V. All rights reserved. PII S0 9 2 7- 7 7 5 7 ( 9 8 ) 0 06 5 3 - 0

306

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

for most complex fluids. There is a significant gap between atomistic or molecular simulation and continuum simulation techniques in the literature. A concentrated emulsion is a typical complex fluid. It occurs in a wide range of industrial applications including food processing, enhanced oil recovery, spray painting, emulsion polymerisation, advanced materials (e.g. composed ceramics) manufacture, and also in nature as biological fluids. One of the fundamental problems is to understand and to predict the rheological properties of concentrated emulsions from their micromechanics. Although Einstein’s work [1] for suspensions of rigid particles was extended by Taylor [2] to dilute drops in the early 1930s, the fundamental understanding of concentrated emulsions under flow is still very limited [3,4]. The difficulties are seemingly insuperable with a purely analytical approach because drops are deformable bodies in flow and there are long-range many-body hydrodynamic interactions between drops. The rheology of concentrated emulsions depends on the volume fraction of the dispersed phase w, the viscosity ratio of the dispersed phase to the continuous phase l, and the dimensionless capillary number. The capillary number is defined as the ratio between the viscous force and the surface tension, Ca=g ka/C, where g is the continuous phase 0 0 viscosity, k is the imposed nominal shear rate, a is the undeformed drop radius, and C is the interfacial tension coefficient. In addition, the concentration of surfactants, the spatial distribution and the diffusion process of surfactant molecules close to interfaces play a very important part in determining the complicated behaviour of drops including deformability, coalescence and disruption. Furthermore, viscoelasticity of two-phase materials can give a complete new dimension to the already complicated problem. An ideal simulation technique should be capable of tracking the highly non-linear and non-uniform interfacial dynamics of drops under complex flows in a unified way. The most popular numerical technique for studying the problem is the boundary-integral method that was pioneered for modelling viscous drops in shear flow by Rallison and Acrivos [5,6 ]. The interfaces of drops are described by a set of markers as Lagrangian fluid points. The velocity

of the markers is obtained by solving a Fredholm integral equation of the second kind. The computational task easily becomes prohibitive as the number of drops increases. For unbounded periodic shear flow, Pozrikidis [7] and Li et al. [8] have improved the method by applying the Ewald summation technique in the evaluation of the Stokeslet for studying the ordered drops with the restriction of l=1. The much more realistic simulations of dense emulsions had not been achieved until very recently when the tabulation procedure was introduced in the Ewald sums independently by Li et al. [9] for two-dimensional (2-D) random drops with l=1 and by Loewenberg and Hinch [10,11] for three-dimensional (3-D) random drops with various l. Stone and Leal [12] using the boundary-integral method in conjunction with the time-dependent convective-diffusion equation for surfactant transport have studied the effects of surfactants on drop deformation and break-up. Using a phenomenological interface constitutive equation, Pozrikidis [13] has incorporated the effects of interfacial viscosity in the boundaryintegral formalism. However, there are some limitations in the boundary-integral approach: (1) it is not applicable to viscoelastic (or non-linear) fluids; (2) it would be difficult to handle the Marangoni effect or to take into account of the spatial inhomogeneous distribution of surfactants; (3) maintenance of boundary elements or grids in a consistent way is very difficult if coalescence and disruption processes must be taken into account. It should be noted that, as the number density of drops increases (interfacial area increases dramatically), the difference of the required number of the elements between the boundary element methods and the finite-element methods becomes smaller. Hence, the significance of the boundary-element method in saving CPU cost is reduced. The development of alternative methods is still very much needed at present. Mo and Sangani [14] have extended the Stokesian dynamics [15] technique to study numerically 3-D random drops under shear in the limiting case of very small capillary number for which the deformation of drops is negligible. Nobari and Tryggvason [16 ] have solved the full Navier–Stokes equations using a front-tracking finite-difference method [17] that

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

allows a fully deformable fluid interface. 3-D simulations of the collisions of two drops have been studied in detail. Our approach to the problem is based on a novel Lagrangian–Eulerian strategy. It combines the Lagrangian nature of the molecular dynamics technique for tracking microstructure evolution of complex fluids under flows and the Eulerian nature of most finite-element (or finite-difference) methods for efficiently solving continuum partial differential equations into a unified numerical scheme. The scheme can be summarised as follows. Given a flow geometry, microstructure of complex fluids is represented by a grid at physical time t (starting from t=0). The solution of time-dependent physical quantities, such as velocity, pressure, stress and strain on the grid, is obtained by solving a set of governing equations in an Eulerian part of the simulation. The grid is then adapted along with moving boundaries in a Lagrangian manner by a given physical time step dt. The iterative scheme is repeated again at a new physical time t∞=t+dt until the flow reaches steady state or becomes fully developed. It is completely free from the above limitations of the boundary-element methods and has been applied in modelling a range of complex fluids such as concentrated suspensions [18], random foams [19] and continuum viscoelastic fluids [20–23]. The diversity of these applications illustrates the generality and the flexibility of the method, and indicates that the method is quite promising for modelling the dynamics of structured fluids. Its advantage for modelling emulsions under flows will be apparent after we illustrate the details of the scheme. Voronoi tessellation has a geometrical property of switching the neighbouring cells’ connection without loss of each cell’s volume (or area in 2-D). We preserve Voronoi cells as Lagrangian fluid elements in the present numerical scheme. The difficulty acknowledged in the previous formulation [20,21] that arises from the mismatching between the number of equations and variables has been eliminated. Our eventual aim is to study the interfacial dynamics of nonlinear multiphase fluids. For the purpose of presenting the scheme and demonstrating its feasibility, we focus in the present paper on 2-D simulations of viscous dense emulsions under

307

shear. The development of 3-D computer code is also under way. The computational model and numerical scheme are illustrated in Section 2. The results are reported in Section 3 and we end with a summary in Section 4.

2. Computational model for emulsions Voronoi tessellation [24] is used as an adaptive grid in our computational model. It is defined as a subdivision of space determined by a finite set of distinct points: each point (centre) has associated with it the region of space nearer to that point than to any other. In the 2-D case, the Voronoi tessellation is a set of polygonal cells in the plane surrounding each centre such that the edges are perpendicular bisectors of the lines between all neighbouring pairs of points. Thus, the vertex of a Voronoi cell is the circumcentres of triangles defined by adjacent triplets of centres, and the centres themselves form the Delaunay triangulation. The advantages of using Voronoi cells as Lagrangian fluid elements are as follows: a) Voronoi polygons map the space defined by the centres and bounding edges uniquely, none of the polygons overlaps; b) the area and edge length of each polygon change continuously when the centres move. The dual lattices, Voronoi centres and vertices, are sketched in Fig. 1(a). When a centre moves inside of the circumcircle of the triangle formed by other centres, both the centres and the vertices must be reconnected as shown in Fig. 1(b) by satisfying the definition of the Voronoi tessellation. Thus whenever the Voronoi grid distorts sufficiently to affect numerical accuracy and convergence, the grid may be automatically and uniquely adapted. 2.1. Emulsions on Voronoi tessellation In principle, we may colour the Voronoi tessellation in order to designate different phases of complex fluids. Each individual cell follows its designated physical laws according to its colour. However, great care must be taken with the interface dynamics in terms of both physics and grid reconnection. For two-phase fluids, the local

308

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

Fig. 1. (a) The Delaunay triangulation (sold line) and the Voronoi tessellation (bold line); (b) Reconnection from AC to BD.

volume fraction y of the dispersed phase at each Voronoi cell is used as an order parameter for designating the continuous phase (y=0) and the dispersed phase (y=1) in our model. The common Voronoi edges of the neighbouring cells, which have different values of the order parameter, represent physical interfaces. Their vertices are Lagrangian points for tracking the microstructure evolution of drops under flow. The other Voronoi edges in bulk phases have no physical significance and are only for the discretisation of continuous phases in the usual sense of finite elements. Notice that the reconnection of vertices by the Voronoi rule provides a natural geometric operation in our scheme for the realisation of coalescence and rupture of drops, though physical criteria need to be determined. In the present paper, our simulations are limited to the regime where the time scale for coalescence and rupture of drops is much longer than the rheological time. The coalescence and rupture of drops are therefore not considered hereafter. For some applications involving interfacial diffusion of materials, a continuous function of the order parameter could be adopted and its local value would be a part of the global solution. The development along this direction will be published elsewhere. Up to now, the problem of dynamic modelling concentrated emulsions under flow has been transformed to the problem of tracking the evolution

of the Voronoi cells and the Lagrangian vertices in a given flow geometry and boundary conditions. At any physical time t, the solution of timedependent physical quantities such as velocity, pressure, stress and strain on the grid, can be obtained in the Eulerian part of the simulation by solving a set of the governing equations given below: V · v=0

(1)

and V · [s−pI]=[df ] (2) interface where v and p are velocity and pressure fields respectively; I is a unit tensor; [df ] is the interface force arising from the interfacial tension between the two phases. We consider Newtonian fluids in the present paper. The stress s is simply given by the Newtonian law, s=g(y) (Vv+Vv†)

(2a)

where g(0)=g and g(1)=lg are the viscosity of 0 0 the continuum phase and the dispersed phase respectively. Thus Eq. (2) can be rewritten as g(y)V2v−Vp=[df ] (2b) interface The generalisation of this method to viscoelastic fluids can be achieved easily as has already been done in Refs. [20,21].

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

2.2. Numerical discretisation on Voronoi tessellation A staggered grid commonly used in computational fluid dynamics for the treatment of the primitive variables is adopted in our numerical scheme. Velocity is designated at the centres of the Voronoi cells as well as at the vertices. Pressure is assumed to be uniform within each Voronoi cell and is defined at the cell centre. The numerical operators corresponding to the gradient and divergence operators are constructed in the following way. Let us focus our attention on the quadrangle, e(ij ) as shown in Fig. 2. Its four nodes are the cell centre i, one of its vertex j, and two other point b and c located at the middle of the edges connected to the vertex j respectively. A physical quantity at the middle points is simply given by the average quantities of the adjacent vertices. Using the Gauss divergence theorem, the average velocity gradient (Vv) over the element e(ij ) is approximated by e(ij) the line integral (in 2-D case) as (Vv)

e(ij)

=

1

P

A e(ij) k

n v dl = k k k A

1

∑ n v Dl k k k e(ij) k (3)

309

where k denotes one of the four edges of the quadrangle e(ij ), i.e. ib, bj, jc and ci as shown in Fig. 2; Dl is an edge length of the edge k; n is an k k outward unit vector normal to the edge k and v k is the edge velocity given by the average velocities of its adjacent nodes; A is the area of the e(ij) quadrangle e(ij ). The average Newtonian stress of the element is readily expressed by s

e(ij)

=g(y)[Vv+Vv†]

e(ij)

(4)

We design Eq. (2) to be satisfied at all Voronoi centres where the uniform pressure within each cell is assumed and there is no surface tension force at the centres. Let c(i ) stand for the Voronoi cell represented by the center i. Since the pressure p is a constant in c(i ), the force balance Eq. (2) can be rewritten as

P

1 [V · s] = (n · s ) dl c(i) A j j j c(i) 1 = ∑ [n +n ] · s Dl =0 bj jc e(ij) e(ij) j A c(i) j

(5a)

where the summation is taken over the quadrangles which share the centre i; n and n are the outward bj jc unit vectors normal to the edges bj and jc respectively and their length Dl =|bj|+|jc|; A is the j c(i) area of Voronoi cell i. It is also required that Eq. (2) is satisfied at all vertices. Let e( j ) stand for the element that is composed of the quadrangles e(ij ) with a common vertex j. Thus the equation can be expressed at the vertex j as [V · s−pI]

e(j)

= =

1 A c(i) 1

P

n · [s −p ] dl i i i i

∑ [n +n ] · [s −p ]Dl ci ib e(ij) e(ij) i i A c(j) i =J( j )[df ] (5b) interface

Fig. 2. The fluid elements used in the present numerical scheme.

where the summation is taken over the quadrangles which share the vertex j; A is the total area of e(j) the quadrangles; n and n are the outward unit ci ib vectors normal to the edges ci and ib respectively and their length Dl =|ci|+|ib|. The function i

310

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

J( j )=1 if the vertex j on interface, otherwise J( j )=0. The surface tension force on the vertex is approximated by

relates the velocities and positions of image points to the velocities and positions of real points in the simulation box as

2 [df ] =C ∑ t interface k k

vimage =vreal −n kL x x y y vimage =vreal y y rimage =rreal −n L −n kL t x x x x y y and

(6)

where the sum is over the tangent unit vectors of the physical edges pointing away from the vertex. The interfacial tension coefficient C is a constant in the present model. The cell pressure p in i Eq. (5b) must be adjusted for satisfying the incompressible condition for every cell c(i ), i.e. 1 [V · v] = ∑ (n +n ) · v Dl =0 c(i) A (i ) bj jc j j j c

(7)

where the summation is same as in Eq. (5a) and v is the velocity of the vertex j. It should be noted j that the old ad hoc formulation of the strain rate tensor [20,21] that arises from attempting to remedy the problem of mismatching between the number of equations and variables in the previous discretisation scheme is completely eliminated. 2.3. Computational scheme Given an area fraction of drops and their sizes, the initial configuration of drops are prepared from a hexagonal lattice in a periodic simulation box. The Lagrangian vertices and physical interfaces are automatically specified in consequence of colouring the hexagonal cells. The interfacial energy of the resulting two-phase fluids with rough interfaces are then relaxed by solving the interfacial tension driven flow, i.e. Eqs. (5a), (5b) and (7) without imposing any external shear flow (k=0). In every time step, the Lagrangian points are moved according to the solution of the velocity field from Eqs. (5a), (5b) and (7). The flow would not cease until smooth and circular drops are realised in the simulation box. Typical results are shown in Figs. 3(a) and 5(a). After the equilibrium state had been realised, a Couette flow with constant shear rate k=dv /dr x y was applied to the simulation box by using the Lees–Edwards boundary condition[25]. It simply

(8a) (8b) (8c)

rimage =rreal −n L (8d) y y y y respectively at physical time t where n and n are x y integers; L and L are dimensions of the simulax y tion box along the flow direction x-axis and the shear gradient direction y-axis. We made no assumption about the velocity profile within the simulation box and the resulting velocity field is part of the solution of the governing equations. The Lees–Edwards boundary condition has no artificial effect on the physics picture of complex fluids under shear as long as the characteristic length of complex fluids is smaller than the size of the simulation box. At any physical time t, Eqs. (5a) and (5b) for a fixed pressure field are solved by the preconditioned conjugate gradient technique [26 ] and the cell pressure is then adjusted in the Eulerian iteration till the incompressible condition, Eq. (7), is satisfied. The tolerances of satisfying Eqs. (5a), (5b) and (7) are less than 10−5 in all the simulations. For the given physical time step dt, both the fluid elements ( Voronoi cells) and the Lagrangian vertices on the interface are moved according to the solution of the velocity field obtained from the Eulerian iteration. The new shape of drops is calculated by fitting the updated positions of the Lagrangian vertices using the cubic splines. Vertex on interface is then relocated at the intersection of the fitted smooth interface and the bisector of the neighbouring cell centres which share the same vertex and have same order parameter y. The reconnection of other non-interfacial cell centres and vertices follows the standard Voronoi rule as shown in Fig. 1. Each drop is then subjected to a small uniform expansion or shrinkage to eliminate any accumulative numerical errors on its area. The procedure is repeated at the next physical time

311

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

(a)

(b)

(c)

(d)

(e)

(f)

Fig. 3. Two closely interacting drops in unbounded periodic shear flow at strains of (a) 0.0; (b) 5.0; (c) 10.5; (d) 11.5; (e) 12.0 and (f ) 13.0. The results obtained from 400 cells grid and other simulation details are given in the text.

until that the flow reaches its steady state. The time step dt was determined by the required macroscopic strain. The simulation results have no significant difference for the chosen maximum strain less than 1%. However, the Voronoi cells reconnection routine would encounter difficulty if the strain were chosen to be greater than 1%. Therefore, 1% was chosen in our simulations. All the simulation results are reported in the dimensionless form by choosing the initial drop radius a and the inverse of the imposing shear rate 1/k as the units of length and time respectively. The average bulk stress may be obtained by ensemble average [27]

taken as

T P C

1

s = ab A

+

A−∑ Am m

A BD U

∂v a −pd +g ab 0 ∂r b

T GP C A P HU 1

∑ A m

×dA+C

dA

∂v ∂v a+ a −pd +lg ab 0 ∂r ∂r Am b b

BD

t t dl (9) a b Cm where A is the total area of the simulation box;

312

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

A is the area occupied by individual drop m and m C is its surface; t is the unit tangential vector m pointing in the counter-clockwise direction around the interface of each drop; a and b are the component of the stress tensor respectively. The first term in the above equation results from the continuousphase fluid. The second term is the contribution due to the presence of the drops and the summation is over all the drops in the simulation box. Thus the effective viscosity and the first normal stress difference can be readily calculated by

g=

s  xy k

(10a)

and

N1= s −s  xx yy respectively.

(10b)

3. Numerical results and discussions As a test case, we performed the dynamic simulation of two closely interacting drops in shear flow. It is a crucial test for any method that deals with concentrated emulsions. We set area fraction w= 0.2 and Ca=0.1. Two sets of grids were used. One grid has 400 cells and other has only 196 cells. Each drop is composed of 37 and 19 cells respectively. The initial positions of drops on the y-axis were same (see Fig. 3(a)). The shapes and the configurations of the drops at the five successive instants during their first collision are shown in Fig. 3(b)–(f ). The pictures demonstrate that deformable drops can glide past each other. The mesh in the thin film region between the interfaces can be automatically adapted by the reconnection rules. Fig. 4(a) shows the drops intend to move away from each other in the y-axis after successive collisions. Fig. 4(b)–(e) demonstrate the correlation between the Taylor deformation D =(L −L )/L +L ) where L and L are t a b a b a b the length of the long and short axes of the drops respectively, the inclination angle of the long axis of the drops with respect to the x-axis, a, and their effective viscosity, first normal stress difference. Within the range of the statistical error, the results

obtained from the different grids have good agreements. The results agree qualitatively with the results of Loewenberg and Hinch [10,11] in their 3-D simulation. As drops approach, a laddershaped grid automatically forms in the thin-film region between drops. The lubrication force in the thin-film region could be divergent, as the distance between particles becomes much less than the radius of particles. It can actually keep particles away from each other. Without adding any explicit lubrication force in our computational model, which includes the naive approximation to the thin-film flow field by the three points velocities crossing the thin-film and the cell pressures, the numerical results indicate that the model can capture the essential physics of interacting drops (having constant surface tension on interfaces) under shear flow in the regime of Ca=O(1). It can resolve the thin-film flow between closely approaching drops and can prevent the drops from overlapping in the case of a capillary number as low as 0.1. It was rather surprising. We suggest that it is very important to keep the mean distance between the Lagrangian vertices on the interface less than the lateral length scale d of the thin-film which is formed during drops collision. Otherwise, there would be insufficient resolution of lubrication force in the gap of drops. A simple scaling arguments gives d≈aCa1/2 [11] for Ca=O(1). Therefore, the coarse grid was used in the following simulations for 0.1≤Ca≤0.6. In the absence of van der Waals attraction, the asymptotic calculation by Yiantsioa and Davis [28] for Ca<1 shows that drop deformation can prevent drop coalescence. Unfortunately a large number of Lagrangian vertices is required to resolve the details of flow fields for preventing drops from overlapping. It would be impractical to implement our method in a very small capillary number regime unless the explicit lubrication forces were introduced as in the early work for modelling hard disks under shear [18]. We have carried out the simulations of concentrated emulsions over a range of the capillary number (0.1≤Ca≤0.6) with a fixed area fraction w=0.4 and l=1. An equilibrium or start configuration is shown in Fig. 5(a). Each drop is composed of 19 cells and has about 30 Lagrangian

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

(a)

(b)

(c)

(d)

313

(e) Fig. 4. Results as a function of strain obtained by the simulation for the two interacting drops as in Fig. 3: (a) separation of drop centres on the y-axis; (b) average drop deformation; (c) average inclination angle; (c) effective viscosity; (d ) first normal stress difference. The numbered points ($) in the figures correspond to the snapshots given in Fig. 3 from (b) to (f ). The triangle and the plus show the results obtained from the 400-cell grid and the 196-cell grid respectively.

314

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

(a)

(b)

(c)

(d)

Fig. 5. Snapshots of concentrated emulsions under shear obtained by the simulation with l=1, w=0.4 and Ca=0.6 at strains of (a) 0.0; (b) 3.99; (c) 7.75 and (d ) 11.54.

vertices on its interface. At t=0, thirty drops are randomly placed in the simulation box which is composed of 1444 cells and 2888 vertices. The previous studies [9,14] suggest that the initial particle distribution is not significant as long as the number of particles is large enough, i.e roughly more than 20 for the concentration of particles is well below that from random close packing. The pressure solver, which is still the most expensive part of simulations, typically takes about 10–20 s of CPU time per Lagrangian time step on a DEC 500/500 MHz workstation. Notice that the present scheme is about 30 times faster than the simple pressure solver in the early work [20,21]. The previous study of 2-D drops under shear

by Li et al. [9] is for relatively low capillary numbers (0.08≤Ca≤0.233). In their simulations, the drops show a convex shape. We carried out simulations with higher capillary numbers. A typical microstructure evolution of drops in the shear flow for Ca=0.6 is presented in Fig. 5(b)–(e). The deformation of drops depends on the surrounding local flow field. The highly inhomogeneous behaviour of drops and the dumbbell shape of drops are clearly observed. Our numerical technique can track the evolution of interfaces in a natural physical way. Highly deformed drops linked to each other by thin-film regions (characterised by the ladder-shaped grid as shown in Fig. 5(b)) form a number of small clusters along the compressional

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

(a)

(b)

(c)

(d)

315

Fig. 6. Results as a function of strain obtained by the simulations with l=1 and w=0.4: (a) drop deformation; (b) inclination angle; (c) effective viscosity and (d) first normal stress difference.

axis of the shear flow. The clusters are subsequently broken apart along the extensional axis, as shown in Fig. 5(c). Fig. 6 shows the transient behaviour of the drop deformation, the orientation, the shear viscosity and the first normal stress difference. The shear stress overshoots are observed in the cases of the large capillary numbers. Steady state results are achieved in strains of roughly 5–10. Thereafter, the physical quantities fluctuate around the mean values of the steady states. This presents clearly the correlation between the microstructure of drops in the flows and the macroscopic rheological properties. Deformable drops can glide past each other with less resistance than in the case of hard suspensions and hence reduce the shear stress

contribution. The elongation of the drops in the flow direction produces the large first normal stress difference. The fluctuation of the macroscopic quantities in the steady states results directly from the evolution of the interfaces, the forming and breaking of the small size clusters of drops under flow. Fig. 7(a) plots the steady state deformation of drops D and the inclination angle a as a function t of capillary number. It shows that, as the capillary number increases, the deformation of drops becomes larger and the elongated drops align more closely along the flow direction. The slight downwardly concave drop deformation curve indicates that there might be an asymptotical value of the

316

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

(a)

(b)

(c)

Fig. 7. Steady-state results as a function of the capillary number: (a) drop deformation and inclination angle; (b) effective viscosity and (c) first normal stress difference.

average D as the capillary number increases and t subsequently there should be no drop break-up under the simulation condition. This is actually expected, since Rayleigh instability of a drop when it is twice as long as its diameter is a 3-D mechanism. By closely examining our animation of the 2-D drops under shear, it is also observed that the drop rotates from the extensional axis to the compressional axis before being significantly deformed. Fig. 7(b) plots the steady state effective viscosity as a function of capillary number and shows the shear-thinning behaviour. It is mainly due to the orientation of the drops in the flow direction. Drops can deform and can glide past each other and produce much less resistance than

in the case of hard suspensions. Notice that the steady-state effective viscosity of hard disks or spherical drops of Ca=0 and l=2 obtained by the early simulations [17] is 3.75 for area fraction w=0.4. Shear rate has much greater effect on reducing the shear viscosity in the very low capillary number regime than in the high capillary number regime. It is consistent with the early observation that the average deformability of 2-D drops would be stagnated as the capillary number increases. Fig. 7(c) plots the steady state first normal stress difference as a function of capillary number. The dependence of the normal stress on shear rate changes from quadratic to linear as the shear rate increases. It is a much more sensitive

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

quantity than the shear stress to the deformation and alignment of drops in high capillary regime. The linear dependence of the first normal stress difference on the shear rate can be qualitatively understood as follows. For a long drop lying along the x-axis in the steady state, the contribution to the first normal stress difference from the interface of the drop can be calculated by (assuming t #1 x and t #0) y

P

dl(t2 −t2 )#2CL #2Ca2/L x y a b Cm since L L #a2 is a constant. The viscous force, a b which tends to stretch the drop, is of the order of lg kL . It must balance with the surface tension 0 b C, thus L #C/(lg k). We therefore conclude b 0

N132lg a2k. Many drops exhibit a concave 0 shape in this regime, but in 3-D, drops might break up at such a high capillary number. C

4. Summary We have applied a general computational technique to study hydrodynamically interacting 2-D concentrated emulsions under Couette flow. The results of the previous studies by Li et al. [9] and Loewenberg et al. [10] have been confirmed. We carried out simulations at higher capillary number. In this regime, the highly inhomogeneous behaviour of drops and the concave shape of the drop can be clearly observed. The fluid exhibits the strong characteristics of viscoelastic fluids: shearthinning of shear stress and non-zero of first normal stress difference. Deformable drops can glide past each other with much less resistance than in the case of hard suspensions and hence reduce shear stress. We have found that the dependence of the normal stress on shear rate changes from quadratic to linear as the shear rate increases. The proposed new discretisation scheme has completely eliminated the old ad hoc formulation of the strain rate tensor [20,21]. We preserve the favourite choice of using Voronoi cells as Lagrangian fluid elements. Using the preconditioned conjugate gradient technique [26 ], the efficiency of the pressure solver in the Eulerian part

317

of iterations has been greatly improved. Our numerical scheme provides a natural platform for introducing physical mechanism (based on microscopic model ) of the coalescence and rupture of drops into mesoscopic emulsion simulations under flows. The Lagrangian–Eulerian approach can track the details of microstructure evolution in a natural physical way. It is especially useful for mesoscopic modelling of complex fluids in which the characteristic time and length scales are usually several orders of magnitude higher than those for simple fluids. The scheme could be readily applied in studying electrohydrodynamics of two-phase fluids, dynamics of drops in viscoelastic media and many complex flow problems. It can also be applied in solving a general two-fluid model [29,30] including the standard hydrodynamic governing equations, thermodynamic force and diffusive equation for composition field in a set of dynamic equations. Once a free energy functional, a constitutive relation of viscoelastic stress and initial condition have been given, computer simulations of the model can be carried out. The details of such developments will be published in forthcoming papers.

Acknowledgments This work was funded by Kao Institute for Fundamental Research, Kao Corporation, Japan. X.F. Yuan wish to thank Prof. M. Doi, Mrs. K. Uchida, Dr. J. Mino and Dr. K. Hara for their hospitality during his stay in Nagoya and to EPSRC ( UK ) for awarding him an advanced fellowship and research grants (GR/L51010 and GR/M06956).

References [1] A. Einstein, in: R. Furth ( Ed.), Investigations on the Theory of Brownian Movement, Dover Publications, New York, 1956. [2] G.I. Taylor, Proc. R. Soc. London, Ser. A 138 (1932) 41. [3] J.M. Rallison, Annu. Rev. Fluid Mech. 16 (1984) 45. [4] H.A. Stone, Annu. Rev. Fluid Mech. 26 (1994) 65. [5] J.M. Rallison, A. Acrivos, J. Fluid Mech. 89 (1978) 191. [6 ] J.M. Rallison, J. Fluid Mech. 89 (1978) 191.

318

X.-F. Yuan, M. Doi / Colloids Surfaces A: Physicochem. Eng. Aspects 144 (1998) 305–318

[7] C. Pozrikidis, J. Fluid Mech. 246 (1993) 301. [8] X. Li, H. Zhou, C. Pozrikidis, J. Fluid Mech. 286 (1995) 379. [9] X. Li, R. Charles, C. Pozrikidis, J. Fluid Mech. 320 (1996) 395. [10] M. Loewenberg, E.J. Hinch, J. Fluid Mech. 321 (1996) 395. [11] M. Loewenberg, E.J. Hinch, J. Fluid Mech. 338 (1997) 299. [12] H.A. Stone, L.G. Leal, J. Fluid Mech. 220 (1990) 161. [13] C. Pozrikidis, J. Non-Newtonian Fluid Mech. 51 (1994) 161. [14] G. Mo, A.S. Sangani, Phys. Fluids 6 (1994) 1637. [15] J.F. Brady, G. Bossis, Annu. Rev. Fluid Mech. 20 (1988) 111. [16 ] M.R.H. Nobari, G. Tryggvason, AIAA J. 34 (1996) 750. [17] S.O. Unverdi, G. Tryggvason, Physica D 60 (1992) 70. [18] X-F. Yuan, R.C. Ball, J. Chem. Phys. 101 (1994) 9016. [19] X-F. Yuan, S.F. Edwards, J. Non-Newtonian Fluid Mech. 60 (1995) 335.

[20] X-F. Yuan, R.C. Ball, S.F. Edwards, J. Non-Newtonian Fluid Mech. 46 (1993) 331. [21] X-F. Yuan, R.C. Ball, S.F. Edwards, J. Non-Newtonian Fluid Mech. 54 (1994) 423. [22] O.G. Harlen, J.M. Rallison, P. Szabo, J. Non-Newtonian Fluid Mech. 60 (1995) 81. [23] P. Szabo, J.M. Rallison, E.J. Hinch, J. Non-Newtonian Fluid Mech. 72 (1995) 73. [24] P.J. Green, R. Sibson, Comput. J. 21 (1978) 168. [25] A.W. Lees, S.F. Edwards, J. Phys. C 5 (1972) 1921. [26 ] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in FORTRAN, 2nd ed., Cambridge University Press, Cambridge, 1992, p. 77. [27] G.K. Batchelor, J. Fluid Mech. 41 (1970) 545. [28] S.G. Yiantsios, R.H. Davis, J. Colloid Interface Sci. 144 (1991) 412. [29] M. Doi, in: A. Onuki, K. Kawasaki ( Eds.), Dynamics and Patterns in Complex Fluids, 1990, p. 100. [30] M. Doi, A. Onuki, J. Phys. II France 2 (1992) 1631.