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

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-Delta similar to 10(3) 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 lambda. Extensive long-time simulations for N = 100-200 and N-Delta = 1000-2000 are performed up to c = 0.55 and drop-to-medium viscosity ratios up to lambda = 5, to calculate the nondimensional emulsion viscosity mu* = Sigma(12)/(mu(e)(gamma)over dot), and the first N-1 = (Sigma(11) - Sigma(22))/(mu(e)\(gamma)over dot\) second N-2 = (Sigma(22) - Sigma(33))/(mu(e)\(gamma)over dot\) normal stress differences, where (gamma)over dot is the shear rate, mu(e) is the matrix viscosity, and Sigma(ij) is the average stress tensor. For c = 0.45 and 0.5, mu* is a strong function of the capillary number Ca = mu(e)\(gamma)over dot\a/sigma (where a is the non-deformed drop radius, and a is the interfacial tension) for Ca much less than 1, so that most of the shear thinning occurs for nearly non-deformed drops. For c = 0.55 and lambda = 1, however, the results suggest phase transition to a partially ordered state at Ca less than or equal to 0.05, and mu* becomes a weaker function of c and Ca; using lambda = 3 delays phase transition to smaller Ca. A positive first normal stress difference, N-1, is a strong function of Ca; the second normal stress difference, N-2, is always negative and is a relatively weak function of Ca. It is found at c = 0.5 that small systems (N similar to 10) fail to predict the correct behaviour of the viscosity and can give particularly large errors for N-1, while larger systems N greater than or equal to O(10(2)) show very good convergence. For N similar to 10(2) and N-Delta similar to 10(3), the present algorithm is two orders of magnitude faster than a standard boundary-integral code, which has made the calculations feasible.