Discrete element modeling of the microstructure of fine particle agglomerates in sheared dilute suspension

Physica A 412 (2014) 66–83

Contents lists available at ScienceDirect

Physica A journal homepage: www.elsevier.com/locate/physa

Discrete element modeling of the microstructure of fine particle agglomerates in sheared dilute suspension A. Kimbonguila Manounou, S. Rémond ∗ Université Lille Nord de France, F-59000 Lille, France LGCgE, Mines-Douai, GCE, F-59508 Douai, France

highlights • • • • •

The initial agglomerates are generated by the Cluster–Cluster Aggregation (CCA) algorithm. Break-up and restructuring of loose fractal agglomerates under shear are studied by the DEM. The fragmentation number is validated from the hydrodynamical and chemical contexts, and the particle size. The packing density of the agglomerates is studied as a function of the fragmentation number. The particle size and the chemical context do not affect the exponents of the power laws.

article

info

Article history: Received 6 January 2014 Received in revised form 16 June 2014 Available online 25 June 2014 Keywords: DEM Agglomerates Fine particle Shear Fractal dimension Microstructure

abstract The fragmentation and restructuring under shear of agglomerates of fine mineral particles are studied with the Distinct Element Method. The model used takes into account contact forces, van der Waals forces, and hydrodynamic forces computed with the free-draining approximation. A loose initial agglomerate is submitted to a constant shear rate until reaching a quasi-stationary state, where the number, size and structure of fragment of agglomerates can be considered as constant. The influence of shear stress and size of particles on the characteristics of agglomerates at equilibrium is studied. Fragmentation is controlled by a non-dimensional number, depending on the radius of the particles, shear rate and maximal adhesion force. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The agglomeration of fine particles in suspensions plays a crucial role in a lot of industrial applications like pharmaceutical, food processing, water treatment, civil engineering. . . . It has in particular a determining influence on the rheological behavior of suspensions. Indeed, agglomeration leads to the formation of large porous composite particles containing part of the free water of the suspension and modifying also the particle size distribution of the suspended particles. It is very important to know the characteristics of these agglomerates so as to apply rheological models allowing us to predict the suspension viscosity [1,2]. Agglomeration of fine particles in a resting fluid is quite well understood [3] and leads to the formation of loose and deformable agglomerates. However, under shear, agglomerates tend to restructure and change into denser composite particles containing water.

∗

Corresponding author at: Université Lille Nord de France, F-59000 Lille, France. E-mail addresses: [email protected] (A. Kimbonguila Manounou), [email protected] (S. Rémond).

http://dx.doi.org/10.1016/j.physa.2014.06.023 0378-4371/© 2014 Elsevier B.V. All rights reserved.

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

67

In most of the previous studies on the fragmentation of the particle agglomerate in a fluid flow, it has been explained that the fragmentation process is dominated by a static or dynamic balance between the cohesive forces (or the cohesive strength of the agglomerate) and the hydrodynamic force (or hydrodynamic stress). This fragmentation process has already been studied by a lot of experimental [4–8] and theoretical and numerical approaches [9–14]. In a shear flow, some researchers pointed out that the relation between hydrodynamic shear stress characterized by the product of the fluid viscosity by shear rate, and the final (stable) agglomerate size obeys the power law [4,15]. The final size of the agglomerate for a given flow condition is largely dependent on the initial agglomerate structure [9,10]. The exponent of the power law is called fracture. Based on the balance theory between driving and resisting forces, some similar dimensionless numbers have already been proposed in the literature to predict the final agglomerate size [9–11,16–18] but have also been used to develop a breakage kernel [19]. It has already been previously observed that when the final agglomerate size is scaled with the dimensionless number, almost all data are aligned on a single straight line [9,10]. The authors emphasize that some data tend to deviate from the master line which they explain by the non-linearity of the behavior of constituent particles in the simulation [9] or do not explain [10]. If all studies agree on the establishment of a final size of agglomerates, a quite large dispersion is obtained in the literature on the values reported for the fracture exponent of the power law, and the explanation for this dispersion does not enjoy consensus. Many researchers have suggested certain key ideas for reflection on this dispersion. Recently, Eggersdorfer et al. [12] attributed the offset of the data to different primary particle sizes, bond strength and level of description of the hydrodynamic force because the inter-particle force as well as the hydrodynamic force change with the primary particle size and their ratio influences the steady state agglomerate size. Harshe and Lattuada [19] emphasize that the fracture exponent depends on the type and magnitude of the inter-particle interactions taken into account. The mechanisms involved in the restructuring have already been explored. In the field of colloidal science, the restructuring of particle agglomerates has been often reported from a viewpoint of the change in space-filling properties [8,20,21]. Generally, the space-filling properties of agglomerates are expressed by the fractal dimension. Harada et al. [18,22] studied by discrete numerical simulations, the restructuring from a viewpoint of the coordination number for the non-fractal agglomerates in a shear flow. The authors have demonstrated that a loose non-fractal agglomerate increases the coordination number and it turns into a dense agglomerate owing to the fluid stress. However, in the case of fractal agglomerates, the internal connectivity may be changed with its fractal dimension and its packing density kept constant. Becker et al. [11] showed that tangential forces significantly affect restructuring behavior and neglecting tangential forces leads to a strong overestimation of the restructuring effects. The transition from the rotational regime to the restructuring regime for the colloidal aggregates exposed to a shear flow has been investigated by Becker et al. [23]. Their results show that the probability that an aggregate being affected by restructuring changes depends on the numbers of primary particles. The kinetics of fragmentation and the evolution of aggregates’ size and fractal structure have been particularly studied. If the prediction and control of the size of the agglomerates in the shear field are required, this information is not sufficient to precisely characterize their microstructure. Indeed, less attention has been paid to the density of agglomerates which characterizes the quantity of water trapped inside composite particles. Despite the role played by the restructuring process on the structural change of the fractal agglomerates under shear, very little attention has been paid to the characterization of a contact network formed by the constituent particles. Also, the influence of the primary particle size has not been explored in detail. These three points have to be specified anyway if one wants to apply existing rheological models to predict the viscosity of suspensions. The aim of this paper is to develop a DEM simulation of agglomerates in sheared dilute suspensions in order to identify the influence of shear rate, physicochemical context and primary particle size on the size, structure and density of agglomerates. The paper is organized as follows: the DEM model is described in Section 2; simulation conditions and numerical methods used for the characterization of agglomerates are presented in Section 3; results and discussion are presented in Section 4 and compared to the literature. 2. DEM model presentation The DEM was initially developed by Cundall and Stack [24]. Since this seminal work, it has considerably developed and it is now widely used for the simulation of granular materials. In our simulations, an initial loose agglomerate of fine spherical mono-sized particles of radius a is immersed into a fluid and submitted to a constant shear rate. 2.1. Equations of motion

mi

Ii

dVi dt

dω i dt

=

N 

Fij + Fi,h

(1)

j=1

=

N  j =1

Tij + Ti,h

(2)

68

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

where mi and Ii are respectively the mass and moment of inertia of the particle Pi , Vi and ωi are respectively its linear and angular velocities, N is the number of particles that are in interaction with particle Pi , Fij and Tij are respectively the interparticle force and torque exerted by particle Pj on particle Pi , and Fi,h and Ti,h are respectively the force and torque exerted by the fluid on particle Pi . 2.2. Inter-particle interactions Two kinds of inter-particle interactions are taken into account in this study: contact interactions Fij,c and van der Waals force Fij,v . The retardation effect of van der Waals interactions is neglected for the sake of simplicity. The agglomeration induced by Brownian motion is also neglected in all the study. Indeed simulations are performed with particles of than diameter larger  or equal to 0.5 µm with large shear stresses (between 1 and 100 Pa) which leads to Peclet numbers Pe = 6π ηγ˙ a3 /kB T much larger than unity (comprised between 3.57 × 102 and 7.15 × 106 ). In addition, the suspensions studied are dilute and the solid volume fraction φ0 is comprised between 7×10−5 and 7×10−1 . The inertial effects of the fluid are neglected in comparison to viscous effects as low Reynolds numbers are obtained in all our simulations (Rep = 4ρs γ˙ a2 /η, where ρs is the particle density and γ˙ is the shear rate) ranging between 10−3 and 3.31. 2.2.1. Contact interactions The contact force models used in our study are similar to those used in Refs. [25,26]. The contact force Tij,c can be decomposed into normal Tij,n and tangential Tij,t components. The normal component of contact force is computed with the nonlinear Hertz model according to Eq. (3) [27]:

 Fij,n = −

  1/2 2  3/2 E¯ R¯ ξn + γn E¯ R¯ ξn ξ˙n nij 3

(3)

where ξn is the overlap between solid particles Pi and Pj , R¯ = ai aj /(ai + aj ), ai and aj represent the radii of particles i and j   respectively, γn is the normal damping coefficient and E¯ = Y / 1 − ν 2 , Y being Young’s modulus and ν Poisson’s ratio. The tangential force exerted on particle Pi by particle Pj is computed according to the Mindlin Deresiewicz theory [28,29] (Eq. (4)):





Fij,t = −µs |Fij,n | 1 −

min |ξ t |, ξt ,max



1−

ξt ,max

 3/2   tij

(4)

where µs is the friction coefficient, ξ t is the cumulated tangential displacement given by Eq. (5), and ξt ,max is the maximal tangential displacement beyond which gross sliding occurs, given by Eq. (6):

ξt =



t

Vij,t dt

(5)

t0

ξt ,max = µs



2−ν 2 (1 − ν)



ξn .

(6)

The tangential force Fij,t exerts a torque on the particle Pi ; moreover, a rolling friction coefficient is also taken into account inducing a resistant moment [30]. The resulting torque Tij exerted by particle Pj on particle Pi is thus given by Eq. (7): Tij = ai nij ∧ Fij,t − µr ai |Fij,n |

ωi |ωi |

(7)

where µr is the rolling friction coefficient. 2.2.2. van der Waals interactions It is important to compute the interaction force or energy between two surfaces to theoretically understand interfacial problems that are frequently encountered in physics, biology, and physical chemistry [31,32]. The classical Derjaguin–Landau–Verwey–Overbeek (DLVO) theory is widely used to describe such interfacial interactions between charged surfaces in liquid media [33,34]. The DLVO theory characterizes the total interaction energy as the combination of van der Waals and electrostatic double-layer interactions. The latter are neglected in this paper as indicated below. In addition, there are two conventional approaches to the calculation of van der Waals energy: Lifschitz’s quantum electrodynamics approach [35] and Hamaker’s microscopic approach [36]. Hamaker’s approach applied to the macroscopic bodies is remarkably simpler and is used here to compute the van der Waals interaction. It assumes that the van der Waals interaction between two bodies of arbitrary geometry is the sum of the interactions between each molecular/atom pair in the system [36]. We use here Hamaker’s constant value reported by Bergström [37] computed according to the Lifschitz theory [35].

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

69

The non-retarded van der Waals force Fij,v exerts an attraction between particles Pi and Pj which is given by Eq. (8) [36,25]: Fij,v = −

AH 6

64a3i a3j δn + ai + aj



×

δn2 + 2ai δn + 2aj δn



2 nij 2  δn2 + 2ai δn + 2aj δn + 4ai aj

(8)

where AH is the non-retarded Hamaker constant of the interacting media, and δn the distance between the two particles’ surfaces. The range of van der Waals forces has been limited to δmax = 100 nm so as to decrease the number of interactions that have to be calculated for each grain. Indeed, the intensity of van der Waals forces decreases rapidly with the separation distance between the two particles, and we have checked in preliminary simulations that larger values of δmax did not change the obtained results. Also, Eq. (8) shows that van der Waals forces lead to a singularity when two particles are in contact (δn = 0) and a maximum value for Fij,v has been fixed for δn ≤ δmin = 0.4 nm. Actually, the particles’ surfaces are not perfectly smooth and they cannot approach up to a minimum distance [12,9,10]. 2.3. Fluid–particle interactions The hydrodynamic force exerted by the fluid on each grain is computed with the Stokes approximation (Eq. (9)): Fi,h = −6π ηai Vi − Vf





(9)

where η is the fluid viscosity, Vi is the linear velocity of particle Pi and Vf is the fluid velocity at the particle position. This approximation is only valid for a particle Reynolds number smaller than unity. Moreover, the free-draining approximation is employed, which considers that the fluid velocity at a given particle position is not perturbed by the presence of surrounding particles. This approximation has already been used in several DEM studies on fragmentation [12,11] and Eggersdorfer et al. [12] have shown that it did not change the scaling laws obtained. The fluid also exerts a torque on particle Pi due to shear, which is computed with Eq. (10) [38] Ti,f = −8π ηa3i ωi − f





(10)

where ωi is the angular velocity of particle Pi , and f = 1/2∇ × Vf the fluid vorticity. 2.4. Gravity and bouncy forces To determine the inter-particle cohesion, the so-called granular Bond number evaluates the ratio of cohesive forces to non-cohesive body or surface forces as given in Eq. (11) Bog =

Fcohesion Fbody/surface

.

(11)

The granular Bond number of Eq. (11) is the most general form as the cohesive force may include a combination of capillary, electrostatic, or van der Waals attractive forces and any relevant body or surface force [39]. In this particular study, the granular Bond number of Eq. (12) is most appropriate: Bog =

AH

  2 16π ρs − ρf ga2 δmin

(12)

where ρs and ρf are the particle and fluid density, respectively, and g and δmin are respectively the gravitational acceleration and the minimum distance particle surfaces can approach, a distance where the attractive forces are maximum. Here, the granular Bond number is simply the ratio between the van der Waals force and the particle weight due to gravity and buoyancy forces. Particles with higher van der Waals forces than their weight (Bog > 1) are considered cohesive while particles with lower attractive forces than their weight (Bog < 1) are considered non-cohesive. In all the simulations investigated in this study, the granular Bond number values are much larger than unity (1.23 × 102 ≤ Bog ≤ 1.23 × 106 ) for particle diameter comprised between 0.5 and 5 µm and Hamaker’s constant values ranging from 10−20 to 10−18 J. Thus, gravity and buoyancy forces can be neglected in our simulations. 3. Simulations’ presentation 3.1. Simulation cell and boundary conditions An initial loose agglomerate is placed in the middle of a simulation cell of dimensions (Lx = Ly = Lz) depending on the size of the agglomerate (Fig. 1). Most of our simulations have been carried out with Lx = Ly = Lz = 150D, where D is the diameter of particles. Initially, all the particles have zero linear and angular velocities. Two impassable planes are used as lateral boundary conditions (in

70

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

Fig. 1. Simulation cell and boundary conditions.

x = 0 and x = Lx) and are given opposite velocities (Vwall,x=0 = −Umax and Vwall,x=Lx = Umax ) in order to apply to the fluid a constant shear rate γ˙ = 2Umax /Lx. Note that all the simulated systems were dilute enough so that no contact occurred during simulation between lateral walls and suspended particles. Periodic boundary conditions are used in the Y and Z directions. 3.2. Initial agglomerate generation The initial system is generated with the Cluster–Cluster Aggregation (CCA) algorithm developed by Botet et al. [40], which allows producing agglomerates composed of Np = 2n particles. Np /2 agglomerates of two particles with a random orientation are first generated; then these agglomerates are merged progressively until obtaining a single agglomerate of 2n particles according to the procedure illustrated in Fig. 2. At each step of the procedure, the radius of gyration of each agglomerate is computed according to Eq. (13): R2g =

Np 2 1  rj − r0 Np j=1

with r0 =

Np 1 

Np j=1

rj

(13)

where Np is the number of primary particles composing the agglomerate, rj their vector position, and r0 is the vector position of the barycenter of the agglomerate. The agglomerates produced with the CCA algorithm are quasi-fractal. The number of primary particles in the agglomerate (Np ) can be expressed depending on their fractal dimension (Df ), radius of gyration (Rg ) and fractal prefactor (kf ) according to the fractal law given by Eq. (14) [41]:

 Np = kf

Rg a

Df

.

(14)

The fractal dimension of agglomerates generated with the CCA algorithm in our study was 1.80 ± 0.02, which is in good agreement with other studies [3,12]. Since the parent CCA agglomerates are obtained by diffusion, their structure is almost loopless with a mean coordination number equal to 2 for all agglomerates studied in this paper. 3.3. Characterization of agglomerates Different characteristics of agglomerates are computed along the simulation: the mean number of primary particles per fragment ⟨Np ⟩, the mean radius of gyration ⟨Rg ⟩, the mean coordination number ⟨z ⟩, and the mean fractal dimension ⟨Df ⟩ (only calculated for agglomerates containing at least 10 particles). Several approaches exist to compute the fractal dimension of a three-dimensional fractal object, which however often lead to different results [42]. In this work, the iterative box counting method (BCM) [43,44] is used to determine the fractal dimension of agglomerates. In this method, the agglomerate is covered with a 3D cubic meshing, the size (ε ) of which will

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

71

Fig. 2. CCA algorithm used for the production of the initial agglomerate.

a

b

c

d

Fig. 3. Illustration of the box counting method to determine the fractal dimension of agglomerates; (a)–(c) discretizing the agglomerate with different box sizes; (d) logarithmic plot of the box number NB (ϵ) versus the reciprocal of box size ϵ .

be iteratively divided by two until it becomes lower than the size of primary particles. At each iteration, the number NB (ε) of cells of the meshing that are cut by the agglomerate is computed (Fig. 3). For a quasi-fractal agglomerate, the relation between mesh size ε and number of cut cells NB (ε) follows a power law given by Eq. (15): NB (ε) ∝

 Df 1

ε

.

(15)

The number of cells that are necessary to completely cover the agglomerate is analyzed as a function of the mesh size. A schematic diagram of the double-logarithmic plot of the reciprocal of box size ϵ and the cell number NB (ε) is presented in Fig. 3(d). The slope of the linear segment of the resulting curve corresponds to the fractal dimension Df of the agglomerate. The floc packing density is defined as the ratio between the total solid volume (Vs = 4/3π a3i Np ) occupied by the primary particles and the total effective volume Vt according to Eq. (16):

γ =

Vs Vt

.

(16)

The total effective volume of agglomerates is very difficult to evaluate accurately. We have therefore considered that agglomerates were approximately spherical, and we have considered two different values for the radius of this sphere: (i) the radius of gyration of the agglomerate Rg according to Eq. (13) (density noted γg ), (ii) the radius of the circumscribed sphere (density noted γc ). Finally, at the end of each simulation, the mean packing density of agglomerates is computed. We evaluate the overall packing density of agglomerates in suspension by calculating a mean packing density weighted by the mass of agglomerates

72

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83 Table 1 Parameters used for the simulations. Micromechanical parameter

Value

Unity

Young’s modulus (Y ) Poisson’s ratio (ν ) Friction coefficient (µs ) Rolling friction coefficient (µr ) Normal damping coefficient (γn )

107 0.29 0.30 0.002 2 × 10−5

N/m2 – – – –

0.01–100 (1.02) 0.4 100

J nm nm

1.0 0.001 5–100 (100)

g/cm3 Pa s Pa

512–4096 (1024) 0.5–5 (0.5) 2.65 10−10

–

Physicochemical parameters Hamaker constant (AH )(×10−20 ) Minimum separation distance (δmin ) Maximum separation distance (δmax ) Hydrodynamics parameters Fluid density (ρf ) Dynamic fluid viscosity (η) Shear stress (τ = ηγ˙ ) Numerical parameter and particles properties Number of particles (Np ) Particle size (d) Particle density (ρs ) Time step, 1t

µm g/cm3 s

⟨γ ⟩ as shown in relation (17):   nf nc , i  1  4 3 1  ⟨γ ⟩ = π ai · m c ,i m f i =1

Vt ,i i=1

(17)

3

where nf is the number of fragment agglomerates in the suspension, nc ,i is the number of primary particles of agglomerate i, and Vt ,i is the total volume of the agglomerate, mc ,i is the mass of agglomerate i and mf is the total mass of agglomerates in the suspension. 3.4. Simulation parameters Table 1 presents the main parameters used for our simulations; default values are indicated in brackets. 4. Results and discussion 4.1. Dynamics of fragmentation 4.1.1. Kinetic analysis Fig. 4 presents the evolution of the mean number of primary particles per agglomerate ⟨Np ⟩, of the average fractal dimension ⟨Df ⟩, of the mean prefactor of the fractal scaling law ⟨kf ⟩ (Eq. (14)) and of the mean coordination number ⟨z ⟩ for an initial agglomerate composed of 1024 particles of diameter d = 0.5 µm. The mean number of primary particles per agglomerate and the mean coordination number are normalized by their initial values (z = 2.5). The ⟨z ⟩, ⟨Df ⟩ and ⟨kf ⟩ values are the average weighted by the mass of the agglomerates. The fractal prefactor kf has been determined according to the alternative expression (Eq. (18)) given by Lee and Kramer [45] which is similar to the general expression previously developed by Gmachowski [46,47] to correlate the fractal dimension and the fractal prefactor of the agglomerates:

 kf =  1.56 −



1.728 −

Df 2

2

  − 0.228

Df + 2 Df

Df 

(18)

where Df is the three-dimensional fractal dimension computed by the box counting method (Section 3.3). The validity of Eq. (18) has already been successfully investigated in the literature using computer simulations and experimental methods [45,47]. The evolution of the initial agglomerate under shear can be divided into four phases. The kinetic evolutions of the average fractal dimension ⟨Df ⟩ and the mean fractal prefactor ⟨kf ⟩ are the same (Fig. 4(b)), only the kinetic of the fractal dimension will be discussed further in the manuscript. Both variables are perfectly correlated by the relation (18). In a first phase of ‘‘elongation’’ (0 ≤ γ˙ t < 1.13), the agglomerate is not fragmented (⟨Np ⟩ = 1024). The mean coordination number remains nearly constant at about 2.5, showing that the agglomerate is essentially composed of linear

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

73

Fig. 4. Evolution of ⟨Np ⟩(t )/Np , ⟨Df ⟩(t ), ⟨kf ⟩(t ) and ⟨z ⟩(t )/z as a function of dimensionless time for an initial agglomerate of 1024 particles of diameter 0.5 µm submitted to a shear stress ηγ˙ = 100 Pa.

‘‘tails’’, almost loopless. The fractal dimension initially equal to 1.79 remains approximately constant, showing that almost no restructuring occurs. The configuration of the agglomerate changes considerably before break-up with the development of the deformation. We observe essentially a rotation and elongation of the initial agglomerate which is progressively stretched in the direction of the flow (along the streamlines). Such a deformation process of the agglomerate is very similar to the experimental observation by Blaser [7]. Fig. 5 shows a view of the microstructure of the first phase described by the kinetics of microstructure evolution. In a second phase of ‘‘fragmentation+restructuring’’ (1.13 ≤ γ˙ t ≤ 5.32), the mean number of primary particles per fragment ⟨Np ⟩ decreases towards a constant value that will not be changed until the end of the simulation, showing that fragmentation is rapidly completed. At the same time, the fractal dimension decreases down to a minimum value (about 1.4 in this simulation) showing that fragments of agglomerates are composed essentially of linear chains of particles which tend to align along the flow direction. Fig. 6 shows the microstructure of the flocs in the second phase of the kinetics of fragmentation where clusters are aligned along the fluid flow.

74

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

Fig. 5. Microstructure of the parent agglomerate in the elongation phase.

Fig. 6. Microstructure of flocs during the phase of fragmentation.

In this second phase, there is a competition between fragmentation and restructuring effects. Indeed, the phenomenon of fragmentation tends to decrease the mean coordination number, because some particles which were initially surrounded by two particles are located after fragmentation at the end of the particles’ chain. On the other hand, some particles are locally rearranging in the agglomerates to find a local more stable position which tends to increase the mean coordination number, this increase clearly demonstrates the restructuring initiation (Fig. 6). The restructuring phenomenon seems to be dominant on changing the mean coordination number of the agglomerates. However, separation of the initial agglomerate in chains that are roughly linear tends to reduce the average fractal dimension despite the fact that locally, the connectivity of some particles increases. For the change of the average fractal dimension of the agglomerates, it seems that the fragmentation effect is predominant, that is, having the massive chains (Fig. 6) which possibly contain some particles with a high coordination number, decreases nevertheless strongly the average fractal dimension of the agglomerates. Our results therefore show that the increase of the mean coordination number of the agglomerates does not necessarily lead to an increase of the fractal dimension in this stage of the restructuring process. Indeed, these two quantities, regarding the microstructural characterization of the agglomerates, do not have the same physical meaning. Both parameters should be used to characterize their internal structure as they provide the complementary information. In a third phase of ‘‘restructuring’’ (5.32 < γ˙ t ≤ 28.50), at this stage, the fragmentation ceases, the mean number of primary particles per agglomerate remains constant, then isolated fragments are restructured and some of them are rounded with the rotation of flow individually. As a result, each fragment gradually loses the fractal appearance their parent aggregate had. The deformation and the formation of new bonds between neighboring particles occur due to the relative motion of constituent particles. The restructuring is continued, and as a result, the packing state changes from loose to dense. Indeed, the average fractal dimension and the mean coordination number increase significantly, demonstrating that the shear changes the structure of the agglomerates. Finally (28.50 < γ˙ t ≤ γ˙ t∞ ), the fractal dimension of agglomerates tends towards a constant value indicating that a quasi-steady state has been reached. The mean coordination number ⟨z ⟩ still increases a little in this zone, but this increase is very slow. This suggests that the state reached at the end of some simulations is not exactly a steady state. However, as the mean number of primary particles and the fractal dimension both have stabilized, and as the increase in ⟨z ⟩ is very slow, we consider in the following that simulations obtained after a shear of 100 dimensionless times are stabilized. Fig. 7 shows

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

75

Fig. 7. Microstructure of flocs at the end of the simulation (steady state).

(a) Mean number of particles per floc, ⟨Np ⟩.

(b) Average fractal dimension, ⟨Df ⟩.

Fig. 8. Evolution of ⟨Np ⟩ and ⟨Df ⟩ as a function of dimensionless time for different initial CCA agglomerates containing 512, 1024, 2048, and 4096 particles of diameter 0.5 µm submitted to a shear stress ηγ˙ = 100 Pa.

the microstructure of the agglomerates at steady state. This physical behavior is in very good agreement with previous studies [12,9,10]. 4.1.2. Influence of the initial agglomerate size The influence of the size of the initial CCA agglomerate on the kinetics of fragmentation has also been studied (Fig. 8). We can observe that the time of the first fragmentation globally increases when the size of the initial agglomerate decreases. The shear gradient exerts a linear force per unit length on the ‘‘chains of particles’’ composing the agglomerate. The tension force between linked particles located somewhere in the middle of the agglomerate will therefore be larger for large agglomerates and will reach the maximal adhesion force sooner which leads to a faster break-up. However, Fig. 8 shows that the size of the initial agglomerate does not change significantly the mean number of particles per agglomerate ⟨Np ⟩ and their average fractal dimension ⟨Df ⟩ once the steady state has been reached. Most of the simulations presented in this paper have therefore been carried out with initial agglomerates composed of 1024 particles. 4.1.3. Influence of shear stress Fig. 9 presents the evolution of the mean number of particles per agglomerate ⟨Np ⟩ and the mean radius of gyration ⟨Rg ⟩ with increasing values of the shear stress. The time-lag for break-up and the size of agglomerates both decrease when the shear rate increases, which is in good agreement with previous studies [4–6,12,9,10,48]. In addition, we show an increase in the mean radius of gyration in the deformation phase (increase of the space occupied by the initial agglomerate, i.e. its initial structure and shape are changed over time due to its stretching in the direction of flow) to a peak at the time of fragmentation before decreasing and reaching an asymptotic value which depends strongly on the hydrodynamic shear stress. Scaling laws relating the mean radius of gyration ⟨Rg ⟩ or the mean number of particles per fragment agglomerates ⟨Np ⟩ to the shear stress will be discussed below and compared to the literature. 4.2. Dimensional analysis We have shown in Section 4.1 that the break-up of agglomerates occurs before their restructuring, when the mean coordination number is still equal to 2. The break-up is therefore related to the competition between the resistant force

76

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

(a) Mean number of particles per floc, ⟨Np ⟩.

(b) Mean radius of gyration, ⟨Rg ⟩.

Fig. 9. Evolution of ⟨Np ⟩ and ⟨Rg ⟩ with dimensionless time for different shear stresses (initial agglomerate of 2048 particles of diameter 0.5 µm).

(maximal adhesion force) and the driving force (due to shear stress) acting on a pair of particles (Pi , Pj ) that are very close from one to another. For two spherical particles of the same radius a and smallest separation distance, the van der Waals forces given by Eq. (8) can be rewritten as Eq. (19): Fa =

AH a 12δn2

.

(19)

The traction force exerted between the two particles comes from the difference between the hydrodynamic forces exerted on particles Pi and Pj (Eq. (20)):

  1Fi,h = −6π ηai Vi − Vf ,i − Vj + Vf ,j ≈ 6π ηγ˙ a2 sin (2θ )

(20)

where Vi ≈ Vj are the velocities of the two particles moving together in the fluid, Vf ,i and Vf ,j are the velocities of the fluid at the position of the particle’s centers, and θ is the angle between the fluid direction and the center to center vector. The driving force maximizes for θ = 45° and the break-up of agglomerate should therefore scale with the dimensionless number developed below. Hence, the fragmentation number (G) is defined as the ratio between the hydrodynamic forces (Eq. (20)) and the maximal physicochemical forces (Eq. (19)) according to Eq. (21):

 G=

2 72π ηaγ˙ δmin

AH



.

(21)

The maximal adhesion force corresponds to the van der Waals forces computed for the minimal separation distance between the two particles, δmin . This relationship is similar to the flow number used by Higashitani and Iimura [9] and Higashitani et al. [10] or the dimensionless number defined by Becker et al. [49]. For a given mono-sized suspension, the maximum attractive force is proportional to the particle’s radius; the fragmentation should in that case only depend on the product aγ˙ . 4.3. Influence of the fragmentation number G at steady state The dimensionless number G (Eq. (21)) has been validated on the specific agglomerated suspensions of mono-sized particles, constructed so as to that it is kept constant. The simulation conditions and results are detailed in the Appendix. Fragmentation is controlled by G, depending on the radius of the particles, shear rate and maximal adhesion force. The latter have shown that the behavior of a suspension of mono-sized particles of radius λa submitted to shear stress τ = ηγ˙ should be the same as the behavior of a suspension of mono-sized particles of radius a submitted to shear stress λτ . This number can now be varied in order to establish scaling laws allowing the determination of the main characteristics of flocs in agglomerated mono-sized suspensions. 4.3.1. Presentation of simulations Table 2 presents the various conditions of simulations. Three cases are defined to study the influence of the fragmentation number G. In the first case, the fragmentation number varies with the velocity gradient, and the physicochemical context and the particle size are both fixed; in the second case, the physicochemical context is fixed and the number of fragmentation G varies with the particle size and the velocity gradient; finally in the third case, the particle size and the velocity gradient are fixed and the number of fragmentation G varies with the physicochemical context via Hamaker’s constant.

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

77

Table 2 Definition of the different cases studied in our simulations. Simulations

Fixed parameters

Variable parameters

Case 1

Particle size (d) Physicochemical context (AH )

Shear stress (τ )

Case 2

Physicochemical context (AH )

Particle size (d) Shear stress (τ )

Case 3

Particle size (d) Shear stress (τ )

Physicochemical context (via AH with δmin fixed)

(a) Mean number of particles per floc.

(b) Average radius of gyration of the floc.

Fig. 10. Evolution of ⟨Np ⟩ and ⟨Rg ⟩ with the fragmentation number G for varying particle sizes, shear stresses and Hamaker’s constant. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

4.3.2. Scaling laws on floc characteristic size Fig. 10 shows the evolution (a) of the mean number of primary particles per floc and (b) the mean radius of gyration at steady state with the fragmentation number G for various conditions. The black, red and blue lines correspond to the three cases defined in Table 2. These two figures show that the average number of primary particles per floc ⟨Np ⟩ and the mean radius of gyration ⟨Rg ⟩ increase when the fragmentation number G decreases. The flocs obtained in the second case where the physicochemical environment is known and the particle size and velocity gradient are varying, are generally smaller regardless of the fragmentation number. We can notice that if the x-axis is scaled using the fragmentation number G, all data for different primary particle sizes in varying physicochemical and hydrodynamic conditions are aligned for the mean number of primary particles per floc and the mean radius of gyration. The dependence of the mean radius of gyration on the fragmentation number can therefore be expressed with the following power law (22):

⟨Rg ⟩ a

= K1 G−m

(22)

where K1 is a constant depending on the physicochemical and hydrodynamical contexts of the suspension, and on the structure of the parent agglomerate. We find a geometrical exponent m equal to 0.34 ± 0.01, determined by linear fitting of the averaged data of all our simulations, whatever the size of the parent agglomerate, diameter of particles, shear stress and physicochemical conditions used. The exponent m characterizes the balance between hydrodynamic forces that tend to separate particles and van der Waals forces which tend to keep the particles in cohesion within the agglomerate. It was also used in previous studies to classify the agglomerates according to their mechanical behavior, that is to say, the fragmentation mechanisms to which they respond to, according to the stress applied to them. According to Snabre and Mills [50], the exponent m mainly depends on the reversibility of the deformation of the agglomerates under the action of external forces. They found values of m equal to 1/3 for rigid agglomerates and to 1/2 for deformable agglomerates. Table 3 summarizes the value found for m in our simulations and compares it to literature results obtained with initial agglomerates of similar fractal dimension and similar shear conditions. Potanin [51] has shown that the exponent m only slightly depends on the fractal dimension of the parent agglomerate. These results confirm the numerical simulations and experimental results reported by Refs. [4,5,12,52] and theoretical results obtained by Wessel and Ball [17] from an analysis of the balance forces for rigid agglomerates in a simple shear. The m value is generally reported in the range of 0.3–0.5 in the literature [53,51]. The exponent m = 0.34 ± 0.01 obtained in our simulations shows that the floc internal structure is rigid due to the contribution of non-central forces, i.e. the tangential forces [51]. This value is in good agreement with previous experimental [4,6,5] and numerical studies [12,52] but in varying physicochemical contexts, and with varying particle sizes and shear stresses. This comparison shows that the model is able to describe the mechanisms of the fragmentation of the agglomerates of usual systems such as those studied by the authors previously quoted. The fact that the exponent m of the scaling law is close also shows that the fragmentation mode of the agglomerates is the same.

78

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

Table 3 Comparison of scaling laws’ exponents with previous studies. Reference

m (–)

p (–)

d (µm)

Df (–)

Type of study

Flow type

Hunter and Frayne [54] Sonntag and Russel [4,5] Brakalov [55] Torres et al. [16,56] Harshe et al. [52] Potanin [51,53] Wessel and Ball [17] Potanin [57] Higashitani and Iimura [9] Eggersdorfer et al. [12] Higashitani et al. [10]

0.37–0.47 0.35 0.58 0.50 0.35 0.23–0.49 0.33 0.38–0.56 – 0.35 –

– 0.87 – – – – – – 0.817–0.896 0.74 0.72

– 0.14 – 0.1 0.81 – – – 0.2–2.0 0.5 0.2

1.8 2.2 1.8 1.8 1.7 1.8 1.8 1.8–2.6 1.8 1.8 1.72

Experimental Experimental Experimental Experimental/Num. Experimental/Num. Theoretical Theoretical Numerical Numerical Numerical Numerical

Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear Simple shear

This work

0.34 ± 0.01

0.71 ± 0.01

0.5–5.0

1.8

Numerical

Simple shear

Table 4 Floc structure at steady state for the hydrodynamic shear stress ranging between 5 and 100 Pa for d = 0.5 µm and AH = 1.02 × 10−20 J. Hydrodynamics shear stress, ηγ˙ (Pa)

5

10

20

50

100

Average fractal dimension, ⟨Df ⟩ (–) Average prefactor, ⟨kf ⟩ (–)

2.03 1.27

2.01 1.23

2.03 1.27

2.0 1.24

1.96 1.15

A similar scaling law was found for the mean number of primary particles per floc with the fragmentation number G according to Eq. (23):

⟨Np ⟩ = K2 G−p

(23)

where K2 is a constant which depends on the particle size, and the physicochemical and hydrodynamical contexts. The geometrical exponent p of the scaling law was computed in a similar way as m. The fragmentation rate of parent agglomerates is more important for increasing values of the particle diameter, hydrodynamic shear stress and decreasing Hamaker values. The independence of p with primary particle sizes in parent agglomerates is consistent with the literature [9,19,58]. We have found p = 0.71 ± 0.01, independently of the hydrodynamical and the physicochemical contexts and also of the primary particle size in parent agglomerates. This is consistent with the geometrical exponent value p = 0.725 obtained by Higashitani et al. [10] who investigated by numerical simulations the fragmentation of CCA parent agglomerates with fractal dimension (df = 1.72). Eggersdorfer et al. [12] used discrete element simulations to investigate the fragmentation of parent agglomerates of mono-sized particles with ramified structure (Df = 1.8) subject to a simple shear flow. They have found an exponent p = 0.74 ± 0.02. Other numerical and experimental results in the literature lead to values of the geometrical exponent of the same order of magnitude as in our simulations [4–6,9,10,18]. Fig. 11 shows some microstructures of flocs at steady state after a shear at 100 Pa. The microstructures below are obtained for decreasing values of the dimensionless number G. For very small values of the fragmentation number, the flocs and even the parent agglomerates become difficult to break due to the cohesion which is becoming important because the particle size and the shear rate are fixed. However, a major restructuring has been observed resulting in the formation of dense structures regarding their fractal dimension. 4.3.3. Average fractal dimension of flocs Table 4 summarizes the average values of the fractal dimension and the structure factor obtained at steady state after shear (5 ≤ ηγ˙ ≤ 100 Pa) of the parent quasi-fractal agglomerates containing 1024 primary particles of diameter 0.5 µm. The simulations are performed with the default value of Hamaker’s constant given in Table 1. The fractal prefactor is close to 1.2 for all our simulations. In the hydrodynamic conditions studied here, the average fractal dimension of flocs at steady state does not vary much; it is estimated at 2.03 ± 0.02, and can be regarded as independent of the applied shear stress. The average value of fractal dimension is consistent with the results presented in Refs. [12,16] although the assumptions used to calculate this average are different. Indeed, the average may apply to certain categories of floc obeying some criteria defined by the authors (weighted average, for example by excluding the small flocs or systems where the steady state is not completely achieved, etc.). approach which consists in using the fractal law (14) establishing the relationship between Using  a more  traditional  ln Np and ln Rg /a and the data of the suspensions shown in Table 4, we can also determine the morphological parameters (Df and kf ) of the agglomerates. The double-logarithmic plot of these data gives a straight line (Fig. 12). The slope and intercept of the resulting regression line respectively correspond to the fractal dimension and the logarithmic fractal

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

79

Fig. 11. Microstructure of flocs at steady state with different values of the fragmentation number G for τ = 100 Pa.

prefactor of the agglomerate. By applying this consideration, this time we get a fractal dimension (⟨Df ⟩ = 2.08 ± 0.16) characteristic of the suspension. It is possible to analytically evaluate the fractal dimension of all the suspensions studied in this paper by combining the scaling laws (22) and (23) previously established. This leads to the relationship (24):

⟨Np ⟩ =

K2 p/m

K1



⟨Rg ⟩ a

p/m

.

(24)

80

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

Fig. 12. Graphical determination of the average fractal dimension of flocs with different hydrodynamic shear stresses. The error bars correspond to the uncertainty on the fit of the data to find the best straight possible.

Fig. 13. Evolution of ⟨γ ⟩g and ⟨γ ⟩c with a dimensionless time for different shear stresses (initial agglomerate of 1024 particles of diameter 0.5 µm).

Comparing Eq. (24) with Eq. (14), the mean fractal prefactor and the average fractal dimension can be obtained with Eq. (25): −p/m

⟨Df ⟩ = p/m and ⟨kf ⟩ = K2 K1

;

(25)

with the values of p = 0.71 ± 0.01 and m = 0.34 ± 0.01 we have found ⟨Df ⟩ ≈ 2.03 ± 0.01. We see that this value of the average fractal dimension calculated by this analytical approach via the geometric exponents p and m respectively of the scaling laws of ⟨Np ⟩ and ⟨Rg ⟩ is consistent with the average fractal dimension (⟨Df ⟩ = 2.03 ± 0.02) calculated by the BCM method. The value of the fractal dimension is consistent with those obtained analytically from the scaling laws and the BCM approach. Aubert and Cannell [59] have found that DLA (Diffusion Limited Aggregation) clusters with a fractal dimension of Df = 1.75 ± 0.05 restructure after a sufficiently long time to give more compact objects with Df = 2.08 ± 0.05. Finally, using the expression Df = p/m (Eq. (25)) with the experimental data of Sonntag and Russel [4] with p = 0.879 and m = 0.35, we calculated Df = 2.51, and they reported a measured value of Df = 2.5. 4.3.4. Scaling law for the evolution of floc density Fig. 13 presents the variation of floc density as a function of the fragmentation number G for different simulations with d = 0.5 µm and AH = 1.02 × 10−20 J and a varying shear gradient. The mean packing density ⟨γ ⟩ of the flocs increases monotonically with the fragmentation number according to the scaling laws (26):

⟨γ ⟩ ∝ Gℓ .

(26)

The two methods (circumscribed sphere and radius of gyration) used for packing density calculation lead to very different values. Indeed, the density computed with the mean radius of gyration of aggregates is much larger than the one computed with the radius of circumscribed spheres. We can clearly see that the density of agglomerates increases with an increase in fragmentation number. This result suggests that large agglomerates can trap a significant quantity of free water which can be released when the shear stress increases. However, the analysis of ⟨Rg ⟩ and ⟨Np ⟩ scaled with the fragmentation number showed that the increase of fragmentation leads to a decrease in floc size (de-agglomeration). Thus, the production of small flocs clearly contributes to increasing the packing density of the flocs:

γ = kf



Rg a

Df −3

.

(27)

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

81

From the definition of the fractal law connecting the solid volume fraction, the fractal dimension of agglomerates and their radius of gyration (Eq. (27)) [51] and the scaling law describing the evolution of the radius of gyration of the agglomerates with the fragmentation number G (Eq. (22)), we can establish the following relationship:

γ ∝ Gm(3−Df ) .

(28)

Scaling laws describing the mean radius of gyration ⟨Rg ⟩ and the mean number of particles per fragment agglomerates ⟨Np ⟩ with the fragmentation number G (respectively Eqs. (22) and (23)) have been studied and compared to the literature. The geometrical exponents m and p taken from these scaling laws were connected to the fractal dimension of the agglomerates by the relationship (25). Combining this later with the law evolution of the packing density of the agglomerates with the fragmentation number (Eq. (26)), we finally obtain the following relationship:

γ ∝



2 72π ηaγ˙ δmin

AH

3m−p (29)

with the exponent ℓ = 3m − p = 0.31 which is independent of the particle size, and the hydrodynamic and the physicochemical contexts in the conditions examined here. The average floc density calculated by the two methods is correctly described by Eq. (29) with the different constant of proportionality (0.05 for the circumscribed sphere method and 0.32 for the radius of gyration method) with p = 0.71 ± 0.01 and m = 0.34 ± 0.01 determined above. 5. Conclusions A DEM model has been used to study the restructuring and fragmentation of loose agglomerates of mono-sized spherical particles submitted to a laminar shear flow. The model takes into account contact forces, van der Waals forces, and hydrodynamic forces based on the Stokes force and the free-draining approximation. An initial loose agglomerate generated with the cluster–cluster aggregation algorithm is placed in a simulation cell with a Couette geometry and submitted to shear stresses of varying intensity. The size, structure and density of the agglomerates after shearing have been investigated as a function of the physicochemical context and primary particle size in the initial agglomerate. Four phases have been identified during the flow: elongation, fragmentation combined with restructuring, restructuring and finally a quasi-steady state where the size and structure of agglomerates remain almost constant. Simulation results show that the hydrodynamic shear stress brings about the rearrangement of primary particles and consequently the interparticle connectivity changes inside the agglomerates. By this restructuring mechanism caused by the fluid and the relative motion of constituent particles, the loose agglomerates become dense regarding their fractal structure and packing density. Fragmentation is controlled by a dimensionless number, similar to those already proposed in the literature and depending on the radius of the particles, shear rate and maximal adhesion force. The latter has shown that the behavior of a suspension of mono-sized particles of radius λa submitted to shear stress τ = ηγ˙ should be the same as the behavior of a suspension of mono-sized particles of radius a submitted to shear stress λτ . This number has been validated by a series of discrete numerical simulations by varying the particle size, the velocity gradient or the physicochemical context. General relationships have been proposed between structural and morphological characteristics of agglomerates (fractal dimension, mean number of particles, mean radius of gyration and packing density) and the dimensionless number G depending on the shear intensity, the size of particles and the maximal adhesion force. The steady state densities of agglomerates have been studied with two different hypotheses (spherical apparent volume of agglomerates with a radius equal to gyration radius or with a radius of circumscribed sphere) as a function of dimensionless number. Results show that the density of agglomerates increases monotonically with the fragmentation number. A power law relation determined from the fractal theory is found for the variation of packing density with fragmentation number which allows describing very satisfactorily the trends obtained numerically. This law allows us to predict the packing density of the agglomerates and to gain access to the water quantity they trap. Moreover, scaling laws between the mean number of particles per agglomerate, the average radius of gyration of agglomerates and the fragmentation number G have been derived and are in good agreement with the literature. The m and p exponents of the power laws are independent of the particle size and physicochemical context. They characterize the internal structure and the fragmentation fashion of the agglomerates, respectively. Finally, the fractal dimension of the agglomerates close to 2.0 was determined by three methods which give close results. It is independent of the primary particle size in the initial agglomerate and the hydrodynamic shear stress for a given physicochemical context. Acknowledgments Discussions with Professor J.N. Roux of Navier laboratory (Marne-la-Vallée, France) and Dr. M.L. Eggersdorfer of Particle Technology Laboratory (Zurich, Switzerland) are gratefully acknowledged.

82

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83

Table A.5 Simulations with varying particle sizes and shear stresses having a constant G value. Simulations

A1

A2

A3

A4

A5

A6

A7

A8

A9

τ (Pa) d (µm)

100 0.5

50 1.0

33.33 1.5

25 2.0

20 2.5

16.66 3.0

12.5 4.0

10 5.0

5 10

(a) Mean number of particles per floc.

(b) Average fractal dimension.

Fig. A.14. Evolution of ⟨Np ⟩ and ⟨Df ⟩ as a function of a dimensionless time γ˙ tG with varying particle sizes and shear stresses having a constant G value.

Appendix Validation of the dimensionless number G on mono-sized suspensions of varying diameter The dimensionless number G suggests that the behavior of a suspension of mono-sized particles of radius λa submitted to shear stress τ = ηγ˙ should be the same as the behavior of a suspension of mono-sized particles of radius a submitted to shear stress λτ . We have examined this point by performing several simulations with a constant value of G (Table A.5). Fig. A.14 presents the evolution of the mean number of primary particles per floc ⟨Np ⟩ and the average fractal dimension of the flocs ⟨Df ⟩ as a function of the dimensionless time γ˙ tG for different values of G. The dynamics of evolution and final characteristics of all the suspensions are almost identical, which confirms the relevance of the dimensionless number G to characterize the behavior of mono-sized suspensions in a laminar flow. For a given fragmentation number G, the kinetics describing the variation of the average number of primary particles per floc ⟨Np ⟩(t ) can be adequately described by Eq. (A.1):

  ⟨N ⟩p (t ) = min Np ,

A

(Gγ˙ t )α

 + n∞

(A.1)

where Np is the number of primary particles in the parent agglomerate before fragmentation, n∞ is the mean number of primary particles at steady state, and α and A are the fitting parameters of this model. For instance the parent agglomerate consisting of 1024 mono-sized primary particles of diameter d = 0.5 µm, sheared at 100 Pa gives n∞ = 27.46, A = 2.12 × 10−6 , and α = 3.59 ± 0.01. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

D. Quemada, Eur. Phys. J. Appl. Phys. 1 (1998) 119–127. H. Usui, K. Kishimoto, H. Suzuki, Chem. Eng. Sci. 56 (9) (2001) 2979–2989. P. Meakin, J. Sol–Gel Sci. Technol. 15 (1999) 97–117. R. Sonntag, W. Russel, J. Colloid Interface Sci. 113 (2) (1986) 399–413. R. Sonntag, W. Russel, J. Colloid Interface Sci. 115 (2) (1987) 390–395. R.C. Sonntag, W.B. Russel, J. Colloid Interface Sci. 115 (2) (1987) 378–389. S. Blaser, J. Colloid Interface Sci. 225 (2) (2000) 273–284. T. Serra, X. Casamitjana, J. Colloid Interface Sci. 206 (2) (1998) 505–511. K. Higashitani, K. Iimura, J. Colloid Interface Sci. 204 (2) (1998) 320–327. K. Higashitani, K. Iimura, H. Sanda, Chem. Eng. Sci. 56 (9) (2001) 2927–2938. V. Becker, H. Briesen, Phys. Rev. 71 (2009) 061404. M. Eggersdorfer, D. Kadau, H. Herrmann, S. Pratsinis, J. Colloid Interface Sci. 342 (2) (2010) 261–268. A. Gastaldi, M. Vanni, J. Colloid Interface Sci. 357 (1) (2011) 18–30. M. Vanni, A. Gastaldi, Langmuir 27 (21) (2011) 12822–12833. S.W. Horwatt, I. Manas-Zloczower, D.L. Feke, Rubber Chem. Technol. 65 (4) (1992) 805–821. F.E. Torres, W.B. Russel, W.R. Schowalter, J. Colloid Interface Sci. 142 (2) (1991) 554–574. R. Wessel, R.C. Ball, Phys. Rev. A 46 (6) (1992) R3008–R3011. S. Harada, R. Tanaka, H. Nogami, M. Sawada, J. Colloid Interface Sci. 301 (1) (2006) 123–129. Y.M. Harshe, M. Lattuada, Langmuir 28 (1) (2012) 283–292. M.Y. Lin, R. Klein, H.M. Lindsay, D.A. Weitz, R.C. Ball, P. Meakin, J. Colloid Interface Sci. 137 (1990) 263.

A. Kimbonguila Manounou, S. Rémond / Physica A 412 (2014) 66–83 [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59]

83

V. Oles, J. Colloid Interface Sci. 154 (1992) 351. S. Harada, R. Tanaka, H. Nogami, M. Sawada, K. Asakura, Colloids Surf. A 302 (1–3) (2007) 396–402. Volker Becker, Heiko Briesen, J. Colloid Interface Sci. 346 (1) (2010) 32–36. P.A. Cundall, O.D.L. Strack, Discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65. R.Y. Yang, R.P. Zou, A.B. Yu, Phys. Rev. E 62 (3) (2000) 3900–3908. S. Rémond, Physica A 389 (21) (2010) 4485–4496. N.V. Brilliantov, F. Spahn, J.-M. Hertzsch, T. Poschel, Physica A: Statist. Theoret. Phys. 231 (4) (1996) 417–424. P.A. Langston, U. Tuzun, D.M. Heyes, Chem. Eng. Sci. 50 (6) (1995) 967–987. R. Mindlin, H. Deresiewicz, J. Appl. Mech. 20 (1953) 327–344. Y.C. Zhou, B.D. Wright, R.Y. Yang, B.H. Xu, A.B. Yu, Physica A 269 (2–4) (1999) 536–553. J. Israelachvili, Intermolecular and Surface Forces, second ed., Academic Press, London, San Diego, 1991. H. Ohshima, K. Furusawa, Electrical Phenomena at Interfaces: Fundamentals, Measurements, and Applications, second ed., M. Dekker, New York, 1998. B. Derjaguin, L. Landau, Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes, Prog. Surf. Sci. 43 (1–4) (1941) 30–59. E.J.W. Verwey, J.T.G. Overbeek, Theory of the Stability of Lyophobic Colloids, Elsevier Publishing Company Inc., New York, Amsterdam, London, Brussels, 1948. E.M. Lifschitz, Theory of molecular attractive forces, Sov. Phys. JETP 2 (1956) 73–83. H. Hamaker, The London–van der Waals attraction between spherical particles, Physica 4 (10) (1937) 1058–1072. L. Bergström, Hamaker constants of inorganic materials, Adv. Colloid Interface Sci. 70 (1997) 125–169. L.D. Landau, E.M. Lifschitz, Lehrbuch der Theoretischen Physik, Verlag Harri Deutsch, 1991. M. Capece, Z. Huang, D. To, Marie Aloia, Charles Muchira, R.N. Davé, A.B. Yu, Powder Technol. 254 (2014) 103–113. R. Botet, R. Jullien, J. Phys. A: Math. Gen. 17 (12) (1984) 2517. B.B. Mandelbrot, The Fractal Geometry of Nature, W.H. Freeman, New York, 1983. J.A. Nelson, R.J. Crookes, S. Simons, J. Phys. D: Appl. Phys. 23 (4) (1990) 465. S. Buczkowski, S. Kyriacos, F. Nekka, L. Cartilier, Pattern Recognit. 31 (4) (1998) 411–418. F. Martinez-López, M. Cabrerizo-Vílchez, R. Hidalgo-Álvarez, Physica A 298 (3–4) (2001) 387–399. Chunwoo Lee, Timothy A. Kramer, Adv. Colloid Interface Sci. 112 (1–3) (2004) 49–57. L. Gmachowski, J. Colloid Interface Sci. 178 (1996) 80. L. Gmachowski, Colloids Surf. A 211 (2002) 197. A. Moussa, M. Soos, J. Sefcik, M. Morbidelli, Langmuir 23 (4) (2007) 1664–1673. V. Becker, E. Schlauch, M. Behr, H. Briesen, J. Colloid Interface Sci. 339 (2) (2009) 362–372. P. Snabre, P. Mills, J. Physique III 6 (12) (1996) 1835–1855. A.A. Potanin, J. Colloid Interface Sci. 157 (2) (1992) 399–410. Y.M. Harshe, M. Lattuada, M. Soos, Langmuir 27 (10) (2011) 5739–5752. A.A. Potanin, J. Colloid Interface Sci. 145 (1) (1991) 140–157. R.J. Hunter, J. Frayne, J. Colloid Interface Sci. 76 (1980) 107. L. Brakalov, Chem. Eng. Sci. 42 (10) (1987) 2373–2383. F.E. Torres, W.B. Russel, W.R. Schowalter, J. Colloid Interface Sci. 145 (1) (1991) 51–73. A.A. Potanin, J. Chem. Phys. 96 (12) (1992) 9191–9200. J.S. Shrimpton, M. Danby, Powder Technol. 228 (2012) 241–249. C. Aubert, D. Cannell, Phys. Rev. Lett. 56 (1986) 738.