Journal of Fluid Mechanics

Shear flow of highly concentrated emulsions of deformable drops by numerical simulations

a1 Department of Chemical Engineering, University of Colorado, Boulder, CO 80309-0424, USA e-mail:;


An efficient algorithm for hydrodynamical interaction of many deformable drops subject to shear flow at small Reynolds numbers with triply periodic boundaries is developed. The algorithm, at each time step, is a hybrid of boundary-integral and economical multipole techniques, and scales practically linearly with the number of drops N in the range N < 1000, for NΔ [similar] 103 boundary elements per drop. A new near-singularity subtraction in the double layer overcomes the divergence of velocity iterations at high drop volume fractions c and substantial viscosity ratio γ. Extensive long-time simulations for N = 100–200 and NΔ = 1000–2000 are performed up to c = 0.55 and drop-to-medium viscosity ratios up to λ = 5, to calculate the non-dimensional emulsion viscosity μ* = Σ12/(μeγ[dot above]), and the first N1 = (Σ11−Σ22)/(μe[mid R:]γ[dot above][mid R:]) and second N2 = (Σ22−Σ33)/(μe[mid R:]γ[dot above][mid R:]) normal stress differences, where γ[dot above] is the shear rate, μe is the matrix viscosity, and Σij is the average stress tensor. For c = 0.45 and 0.5, μ* is a strong function of the capillary number Ca = μe[mid R:]γ[dot above][mid R:]a/σ (where a is the non-deformed drop radius, and σ is the interfacial tension) for Ca [double less-than sign] 1, so that most of the shear thinning occurs for nearly non-deformed drops. For c = 0.55 and λ = 1, however, the results suggest phase transition to a partially ordered state at Ca [less-than-or-eq, slant] 0.05, and μ* becomes a weaker function of c and Ca; using λ = 3 delays phase transition to smaller Ca. A positive first normal stress difference, N1, is a strong function of Ca; the second normal stress difference, N2, is always negative and is a relatively weak function of Ca. It is found at c = 0.5 that small systems (N [similar] 10) fail to predict the correct behaviour of the viscosity and can give particularly large errors for N1, while larger systems N [gt-or-equal, slanted] O(102)show very good convergence. For N [similar] 102 and NΔ [similar] 103, the present algorithm is two orders of magnitude faster than a standard boundary-integral code, which has made the calculations feasible.

(Received April 23 2001)
(Revised August 16 2001)

Related Content