Article pubs.acs.org/Langmuir
Utilizing the Discrete Element Method for the Modeling of Viscosity in Concentrated Suspensions Martin Kroupa, Michal Vonka, Miroslav Soos, and Juraj Kosek* Department of Chemical Engineering, University of Chemistry and Technology Prague, Technicka 5, 16628 Prague 6, Czech Republic S Supporting Information *
ABSTRACT: The rheological behavior of concentrated suspensions is a complicated problem because it originates in the collective motion of particles and their interaction with the surrounding fluid. For this reason, it is difficult to accurately model the effect of various system parameters on the viscosity even for highly simplified systems. We model the viscosity of a hard-sphere suspension subjected to high shear rates using the dynamic discrete element method (DEM) in three spatial dimensions. The contact interaction between particles was described by the Hertz model of elastic spheres (soft-sphere model), and the interaction of particles with flow was accounted for by the two-way coupling approach. The hydrodynamic interaction between particles was described by the lubrication theory accounting for the slip on particle surfaces. The viscosity in a simple-shear model was evaluated from the force balance on the wall. The obtained results are in close agreement with literature data for systems with hard spheres. Namely, the viscosity is shown to be independent of shear rate and primary particle size for monodisperse suspensions. In accordance with theory and experimental data, the viscosity grows rapidly with particle volume fraction. We show that this rheological behavior is predominantly caused by the lubrication forces. A novel approach based on the slip of water on a particle surface was developed to overcome the divergent behavior of lubrication forces. This approach was qualitatively validated with literature data from AFM measurements using a colloidal probe. The model presented in this work represents a new, robust, and versatile approach to the modeling of viscosity in suspensions with the possibility to include various interaction models and study their effect on viscosity.
■
predicts a linear dependence of η on ϕ. This relation breaks down for ϕ larger than approximately 5%. The reason for this is that below this threshold the movement of each particle can be treated individually and then added to obtain the behavior of the whole suspension. For larger values of particle volume fraction, the particles inevitably start to influence each other, and the dependence of η on ϕ becomes nonlinear. In this region, phenomenological equations are widely used, such as that of Krieger and Dougherty9
INTRODUCTION Suspensions of particles in a liquid are common in a large variety of applications such as food, blood, and polymer latexes to name only a few. An essential property of these systems is their viscosity because it largely affects their behavior. Conversely, the viscosity of a suspension is generally affected by its properties. The rheological behavior can be very complex, and for this reason, simplified systems are often preferred because they reduce the number of influencing parameters. The simplest system in this context is the hard-sphere model.1 However, even for this simple system the rheological behavior for high particle volume fractions ϕ is still not completely understood. Experimental measurements of the viscosity in hard-sphere suspensions are often complicated by various effects including the deviation from a strictly monodisperse suspension, which introduces a strong effect on the viscosity.2 Furthermore, the presence of interparticle forces (e.g., the van der Waals attraction) introduces additional complexity into the system.1,2 Therefore, benchmark studies of hard-sphere rheology usually used dispersions of silica particles or measurements of chargestabilized latexes mapping the hard-sphere behavior.3−7 For low values of the particle volume fraction, the viscosity η of a hard-sphere suspension follows Einstein’s equation,8 which © XXXX American Chemical Society
−2.5ϕm ⎛ ϕ⎞ ⎟⎟ η = ηf ⎜⎜1 − ϕm ⎠ ⎝
(1)
or Maron and Pierce10 −2 ⎛ ϕ⎞ ⎟⎟ η = ηf ⎜⎜1 − ϕm ⎠ ⎝
(2)
where ηf is the fluid dynamic viscosity and ϕm is the maximum packing fraction. Received: June 23, 2016 Revised: July 29, 2016
A
DOI: 10.1021/acs.langmuir.6b02335 Langmuir XXXX, XXX, XXX−XXX
Article
Langmuir
Each element is characterized by its mass (mi), position (xi), velocity (vi), and rotation rate (Ωi). The governing equation for the translational motion is Newton’s second law
A more recent theory was derived by Mendoza and Santamaria-Holek using the effective volume fraction11 −2.5 ⎛ ϕ ⎞ η = ηf ⎜1 − ⎟ 1 − cϕ ⎠ ⎝
where c =
1 − ϕm ϕm
d2x i
(3)
dt
.
2
=
Fi mi
(4)
where Fi represents the sum of all forces acting on the discrete element i. The temporal change in the rotation rate (Ωi) is expressed as follows
A large amount of effort has been invested in understanding multibody dynamics at large values of particle volume fractions. Models accurately describing the suspension viscosity beyond Einstein’s equation were developed to the square of volume fraction12−14 and recently even to the cubic power of the volume fraction.15−17 These expressions provide correct predictions for particle volume fractions of up to ϕ ≈ 0.3.15 Other approaches to the rigorous modeling of the rheological behavior of concentrated suspensions usually employ Stokesian dynamics (SD) because this method is well suited to the description of the particle−fluid interaction.18−20 Various studies have shown that SD is able to successfully reproduce the dependence of η on ϕ in agreement with experimental data.18,19 These studies, however, either evaluated the suspension viscosity for a fixed configuration of particles (although random) and did not consider dynamics18,20 or used a short-range repulsive potential to prevent particles from overlapping in a dynamic simulation.19 There is a serious limitation of the methods that have been used for the determination of viscosity in concentrated suspensions. It is known that at high particle volume fractions the viscosity is predominantly determined by lubrication forces.1,18 However, the commonly used two-particle expression for lubrication forces diverges to infinity for particles in close contact. This issue has a serious effect on viscosity, but it is not encountered in SD simulations for the reasons described in the previous paragraph. However, considering systems where particles coagulate and form clusters as a result of their adhesive properties, there must be a mechanism that would allow particles to come into solid body contact. One possibility to explain this phenomenon can be the presence of slip on the surfaces of particles. In our model based on the discrete element method (DEM),21,22 we introduce the correction for slip into the two-particle formulation of lubrication forces.23 Furthermore, we constrain the lubrication forces at a certain distance in order to overcome their divergent behavior. A similar procedure was recently used in a different simulation method.24 This allows us to determine the suspension viscosity in a dynamic simulation using a novel method based on the force balance in a simple shear model. We show that with a reasonable value of the slip length our model is able to reproduce both theoretical predictions (eqs 1 and 2) and experimental data for a hard-sphere system.3,5 Moreover, we demonstrate that the viscosity of concentrated suspensions is primarily determined by lubrication forces because they are the most important source of dissipation in the system.
dΩ i M = i dt Ii
(5)
where Mi is the sum of all torques acting on particle i and Ii is the particle moment of inertia. For a homogeneous solid sphere 2 we have Ii = 5 miR p2 (particle moment of inertia around its 4
center) and mi = 3 πR p3ρp. Hydrodynamics. Forces that arise from the interaction between fluid and particles immersed in this fluid are crucial to the rheological behavior of a suspension. The large applied shear rates cause the Peclet number to be very high, and the effect of the Brownian motion thus can be neglected. Because particles considered in this work are small (
