Mikael
Mohtaschemi
*a,
Antti
Puisto
a,
Xavier
Illa
b and
Mikko J.
Alava
a
aDepartment of Applied Physics, Aalto University, P.O.Box 11100, AALTO, FI-00076, Espoo, Finland. E-mail: Mikael.Mohtaschemi@aalto.fi; Fax: +358 9 855 4019; Tel: +358 50 430 2473
bDepartment d'Estructura i Constituents de la Matèria, Facultat de Fìsica, Universitat de Barcelona, Martí Franqués 1, E-08028 Barcelona, Spain
First published on 31st January 2014
We study a colloidal model based on population balances in the context of complex fluid rheology. Two typical particle microstructure kinetics, orthokinetic, collisions due to shear, and perikinetic, collisions due to Brownian motion, are found to appear at continuum as different flow behaviors – those having monotonic and non-monotonic flow curves, respectively. Solving the colloidal model together with the 1D Stokes equation for laminar, incompressible flow with Couette boundary conditions, allows bridging the gap between the rheological experiments and the microstructural modeling. The analysis of such a model reveals that orthokinetic particle suspensions have a uniquely defined, continuous steady state shear profile, whereas suspensions in which also perikinetic collisions are present, the steady state can be shear banded and non-unique. Thus, the shear banded configurations at a steady state are found to depend on the initial conditions and the collision kinetics of the system. At high shear rates all the studied cases show continuous shear profiles.
Unstable colloidal dispersions are one of the great examples of fluids exhibiting many of the described phenomena.3 These materials are vast in technological applications. Often aggregation is intentionally induced due to the requirements of the process. For instance, water removal through filtration from a small particle containing dispersion can be quite cumbersome. Inducing colloidal aggregation makes the water removal process faster by allowing larger pore size for the filters.4 Furthermore, since aggregation processes are usually reversible,5,6 the colloids can be re-dispersed after the required processing.
The complex rheology related to the change of the internal state in the suspensions is still rather poorly understood. Often, the complex fluid rheology dynamics is studied using phenomenological models.7
These models are, for the sake of simplicity, based on abstract variables without physical counterparts, such as the degree of structuration or fluidity, and thus the relationship between microstructural dynamics and rheology is not explicitly considered. For many applications the internal phase structure is as interesting as the fluids rheological properties.8–11 Here, the relevant models become important. The main issues are the relationship between the microstructure and rheological observables, which could shed light on phenomena such as thixotropy and yield stress which still pose many open questions.12 These are the main goals of the present paper.
In unstable colloidal suspensions modeling the microstructure i.e. aggregation dynamics using population balances has a long standing history.13 The popularity of such models rests on the fact, that, unlike most of their particle based counterparts, they are not computationally extremely expensive, allowing to match their time and length scales to experiments for direct comparison. For instance Runkana et al. studied the aggregation through the perikinetic process due to charge neutralization and polymer bridging aggregation.14,15 Such scenarios have been shown to agree well with some experiments, provided the collision frequency is corrected appropriately using the DLVO theory. Perikinetic aggregation and the related rheology under shear have been studied by many authors.5,16,17 The perikinetic processes are further complicated by the fact that the aggregation process itself is reversible and the size distribution is dynamically formed in a competition between aggregation and aggregate breakage due to the hydrodynamic stresses induced by shear.5,18 In general, particles under shear flow can experience both, peri- and orthokinetic aggregation. Of these, one or the other dominates, depending on the flow intensity and particle properties. To reflect this in a model, the simplest approximation that can be taken is to apply a superposition of the kernels describing both processes.19,20
Apart from the present microstructural model, other common models with a proper physical interpretation used to model complex fluid rheology are the Shear Transformation Zone (STZ)21 and the Soft Glassy Rheology (SGR)22 models. The STZ model describes an amorphous solid system of interacting particles near the jamming concentration. There the key observation is that the shear deformation localizes to certain hotspots. The birth and annihilation of these hotspots, called shear transformation zones, are then modeled at continuum in a bistate model. The SGR model, on the other hand, aims to describe soft glassy materials formed of repulsively interacting soft particles well beyond the random packing volume fraction. The particles traverse through a corrugated potential energy surface by jumping from one potential well to another by external activation. The fluid is modeled via mesoscopic fluid elements each having a certain trap energy. Thus, both of these microstructural models treat highly concentrated disordered systems, whereas here, we deal with the dynamics of aggregating colloidal suspensions, in which the particle concentration is well below the jamming concentration.
In this paper, we discuss the rheology of complex fluids using a Population Balance Equation (PBE) based rheological model. We demonstrate that using different kernels for the aggregation and varying the particle concentration, one can access the rheological characteristics of shear thinning (vanishing yield stress), simple yield stress, and thixotropic fluids. We implement the model in a realistic measuring geometry, and compute the fluid's rheological characteristics as they would appear in experimental rheometry. Furthermore, we discuss the observed complex phenomenology comparing them against recent experimental data obtained with traditional rheometry equipped with local velocimetry.23,24 These experiments are usually performed in the concentric cylinder (Couette) geometry, which is also the one applied in the present analysis. In addition we show how changes in the microstructure map to macroscopically observable effects in particular (transient) shear banding and the associated fluidization process.
(1) |
(2) |
ccutoff = cp(1 − tanh((η − ηc)/ηc). | (3) |
The fragmentation kernel describes the splitting of the aggregates under shear. Normally it is assumed that breakage is mainly due to the hydrodynamic stress exerted on the aggregates due to fluid flow.20,30 In principle, depending on the interaction potential of the aggregates, it could be possible that also Brownian motion would induce some fragmentation. However, this effect would be limited to low shear rates. For the fragmentation kernel, no definite first principles expression is presently known. However, several formulations exist in the literature.31,32 Thus, we apply here a simple formulation, demonstrated to capture, together with the aggregation kernel, experimentally observed behavior of colloidal aggregation18
(4) |
The suspension rheology is assumed to be tied to the aggregate size distribution via the jammed volume fraction,33 which consists of the solid particle volume, and the liquid volume embedded inside the fractally growing aggregates. The suspension viscosity is computed using a constitutive equation, which reads18,34
(5) |
(6) |
For the sake of a realistic starting point we use parameters that were obtained earlier18 by fitting the orthokinetic model against rheological data for latex suspensions.16 In particular κ = 1.82, ϕmax = 0.60 in eqn (5) and dF = 2.3, Fc = 0.3 nN, ρ1 = 1.0 μm. Unless explicitly stated the same values are used throughout the text.
To avoid solving exhaustive number of equations implied by eqn (1) we sample the size distribution using the fixed pivot method by Kumar and Ramkrishna.35 The number of classes used in the population balance equation affects the accuracy of the obtained size distribution and naturally also the computational effort needed. This is especially important for distributions with ϕ close to ϕmax (see eqn (5)). We tested different number and spacing of size classes. These cause small systematic shifts in the stress values. However, other quantities such as the local flow profiles in the Couette are not influenced by the discretization. As a trade-off between accuracy and computational efficiency we employ a spacing of up to 65 classes (10 classes with the number of monomers twice that of the previous size class and 5–55 classes with 10 ρ1 steps) sufficient for describing the jammed states where the largest aggregates appear. The numerical solution of eqn (1)–(6) is conducted using the Adams Moulton method as implemented in the SUNDIALS routines.36 In the following sections the model as explained so far is called the homogeneous shear rate model.
With the homogeneous shear rate model it is possible to simulate different types of suspensions of particles with attractive interaction potentials.14,15 Next, the focus is turned to modeling the behavior of such systems in a rheometer, where the shear has spatial dependence. The following discussion is restricted to the commonly employed concentric cylinder Couette device. In order to do so, we introduce a spatial dimension. By assuming axial and translational symmetry and neglecting end effects, one can reduce the problem to 1D (in radial direction). Solving the Stokes equation for this case, assuming laminar flow and no-slip boundary conditions, results in a function for the local shear rate37
(7) |
When a flow curve of a fluid is experimentally determined, stress and shear rate pairs are obtained. In a Couette geometry these quantities are not homogeneous over the gap, and it necessitates to concentrate on a single position. Usually this position is at the rotor, where the torque, and therefore the stress, can be measured. However, even though the cylinder's angular velocities are known, owing to the fact that the viscosity distribution in the gap is unknown, the shear rate has to be approximated. The usual approximation is to assume a homogeneous viscosity distribution over the gap,38,39 which gives at the inner cylinder
(8) |
(9) |
(10) |
Fig. 1 Steady state flow curves obtained with the model at different concentrations and varying amount of perikinetic aggregation. |
In the steady state flow curves the stresses, including the yield stress, drop with decreasing concentration, following experimental findings of the behavior of unstable colloids.16 In the model, the critical shear rate c is a power-law function of both the concentration and the weight of the perikinetic aggregation, cp. The concentration dependence is observed also experimentally.42 The dependence on perikinetic aggregation is experimentally realized by changing the particle interaction potential.43
The cause of the deviations from the intrinsic flow curve is the heterogeneous shear profile, which is due to the stress variation in the gap induced by the circular Couette boundary conditions; due to the constant torque over the gap, the stress is higher close to the cylinder with a smaller radius. In a shear thinning fluid, the higher stress implies a smaller viscosity, which according to eqn (7) means a higher local shear rate, setting up a self-feeding loop.46 In the orthokinetic case, this shows up only at the small shear rates near the stress plateau, where the shear localization appears.47 However, in the superimposed case a more profound effect is seen appearing over a wider range of shear rates. Depending on the initial state, a static shear band of different width appears in the shear rate profiles. When the fluid is initialized in a state having smaller viscosity than the steady state one, the static shear band has the width (hSB) determined by c in accordance to the simple lever rule48
(11) |
If the initial state has a viscosity above the steady state, the high shear band at the steady state is narrower than that given by the lever rule. This is caused by the interplay between the fluid's transient dynamics and the stress distribution of the device. Thus, also the stress measured at the rotor will have a different value, resulting in the observed hysteresis in the flow curves. Steady state radial shear rate profiles for these cases demonstrating the effect are shown in Fig. 3.
The plots of the local shear rate–stress pairs appearing in the Couette gap reveal further, that in the cases applying the superimposed aggregation kernel, the steady state has four scenarios depending on the measurement protocol as shown in Fig. 4 and 5. The first scenario is the case where the fluid is initially in a fluidized state at high shear rate (exhibiting no shear bands). The local shear rate–stress pairs inside the Couette gap for such cases are shown in Fig. 4(a) under shear rate control and in Fig. 5(a) under stress control. In those cases, indeed, the lowest shear rate of the high shear rate band is exactly the critical shear rate of the intrinsic flow curve. The second scenario is shear rate controlled and starts from a jammed state shown in Fig. 4(b). As mentioned before, in that case the resulting shear rate profile depends on the initial state and on the dynamics of the underlying structure, eqn (1). The observed, this time apparent critical shear rate, is larger than the intrinsic one. Comparison of Fig. 4(a) and (b) shows that the simulation initialized in a fluidized state gives homogeneous shear profiles at smaller shear rates than the ones initialized in jammed states i.e. the apparent critical shear rate is always larger than the real one in agreement with recent observations of rheological hysteresis.49 The third scenario is the one where the Couette is operated in a stress controlled mode and the simulation is initialized in a jammed configuration. For such situations the local shear rate–stress pairs are plotted in Fig. 5(b). Now, the left, low stress branch of the flow curve is favored. The applicability of stress controlled operation in experiments, in contrast to simulations, is limited by the long time required to reach the steady state compared to the corresponding shear rate controlled case. In the fourth scenario, the shear rates are well above the critical shear rate and no shear bands appear regardless of the initial condition and the operation mode.
Fig. 5 Local stress-shear rate values at steady state with the superimposed aggregation kernel at different applied stresses. When initial state is (a) fluidized 0 = 100 s−1 (b) jammed 0 = 10−5 s−1. |
Shear rate transients from jammed states to fluidized ones are shown for the orthokinetic case in Fig. 7 and for the superimposed kernel in 8. In both cases, the flow first localizes at the rotor after which a transient shear band forms, which evolves until the high shear band reaches the stator or develops into a static shear band. In the orthokinetic case, the fluidization occurs through a transient shear banding phase, finally ending in a linear velocity profile at high shear rates, and shear localization at small shear rates. The time-evolution of the stress is shown in the upper panels of Fig. 7 and 8. In both cases, at first the stress decays rapidly, but once the transient shear band is formed a slower relaxation time scale related to transient shear banding appears. Finally, the stress settles, as the transient shear band vanishes or develops into a static one. Apart from the missing initial elastic stress-growth, the stress relaxation patterns show remarkable resemblance to the stress pattern recorded in Carbopol gels.40
The stress controlled simulations starting from the fluidized state follow the same pattern as the shear controlled ones; the structure here evolves homogeneously until the steady state is reached. Starting from the jammed state both orthokinetic and superimposed cases result in transient shear banding similar to the shear controlled simulations. A representative plot of this is shown in Fig. 9 for an imposed stress of 2.6 Pa. Similar to the shear rate controlled simulations, the transient shear bands vanish in the orthokinetic case at all imposed stresses. In the perikinetic case the transient shear bands turn into static ones at imposed stresses close to the yield stress. Note, that in the stress controlled mode, the transient shear bands typically have sharper edges, but otherwise the behavior is very similar to the ones under the shear rate control.
How long does one have to wait in order to obtain an (apparent) steady state solution? In particular this question is relevant for the case of yield stress and heavily shear thinning fluids where the structural changes between the flow and no flow states are significant. Divoux et al. addressed the issue of fluidization experimentally for the case of a simple yield stress fluid (i.e. not thixotropic) using extensive rheometry combined with local velocimetry measurements.40,41 The main conclusions drawn are that the fluidization time (reaching a homogeneous flow profile) can take a considerably long time. Moreover, they found that the fluidization time follows power law scalings, with different exponents for the shear rate and stress controlled modes. In the present model, the fluidization time, here defined as the vanishing of a low shear rate band, follows a similar power-law with both shear rate and stress consistent with some recent approaches.46,50,51
In Fig. 10 the fluidization time for shear rate controlled runs as a function of the engineering shear rate is shown. In the orthokinetic case, the fluidization time can be well fitted with a power law having the exponent α ≈ 1. Changing the initial configuration (to more or less jammed) affects only the pre-factor but not the power-law exponent.
Fig. 10 Fluidization times for shear rate controlled runs with superimposed and orthokinetic aggregation kernels. The dashed lines correspond to power-law fits τf ∼ a−α at high shear rate region. |
The results with the superimposed model complicate this simple picture somewhat. First it should be noted that for applied shear rates where the steady state configuration has a static shear band, the fluidization time is not well defined (of course a timescale to reach a steady state can be defined, however). For intermediate shear rates, for which the steady state is not shear banded, the fluidization time follows a power-law with an exponent that increases with cp (see Fig. 10). At very high shear rates ( ≫ 100 s−1) the orthokinetic and superimposed cases coincide.
Next the focus is turned to the fluidization process with stress controlled runs. The fluidization times as a function of the reduced stress are shown in Fig. 11. Here, following the definitions used in the experiments,40 the reduced stress means simply subtracting the yield stress from the imposed stress values. This way, the fluidization times can be fitted by a power-law with the exponent β ≈ 1.6. However, the exponent is quite sensitive to the used yield stress value (here σc = 1.519 Pa obtained from the intrinsic flow curve has been used). At very high stresses the fluidization time scales with an exponent β ≈ 1. The orthokinetic and the superimposed cases are mostly the same. The difference between these cases is that for the superimposed case there is a range of stresses for which the steady state includes a shear band and thus no fluidization time can be defined. Changing the initial configuration (to more or less jammed) does merely affect the pre-factor but not the exponent. Also the scaling at high stresses with an exponent β ≈ 1 and the scaling at the low stress limit with β ≈ 1.6 are independent of the concentration as in the shear rate controlled simulations.
Fig. 11 Fluidization times for stress controlled runs with orthokinetic and superimposed aggregation kernels. σc = 1.519 Pa has been used. The dashed lines correspond to power-law fits τf ∼ a−β. |
To conclude the fluidization times for the shear rate and stress controlled cases follow power law like behaviour. The yield stress can be obtained from the stress controlled fluidization times. Divoux et al. observed that the ratio of the two fluidization exponents gives the exponent of the Herschel–Bulkley model (σ = σo + ηα) fit to the measured materials (apparent) flow curve.41 In the current model the value of the Herschel–Bulkley exponent scales with the κ-exponent in the viscosity function eqn (5). However, the fluidization time exponents were insensitive to changes in κ. Thus, it is possible to tune the parameters of the system to satisfy the aforementioned relationship. Since the relationship does not hold for arbitrary model parameters, the general validity of this scaling relationship is not accessible with the current model.
The aggregate size distribution inside the Couette gap is shown in Fig. 12 for a shear rate controlled run started from the jammed state. The figure shows the explicit relationship of the microstructure and the rheology during the transient. The system is initially in a jammed state. First, the large aggregates are broken down close to the rotor which decreases the local jammed volume fraction and leads thus to a decrease in the viscosity. This mechanism of positive feedback leads to the formation of a transient shear band, which advances towards the stator. However, the shear band does not completely vanish, but evolves into a static one due to the underlying negative slope of the intrinsic flow curve. The steady state has two clearly different aggregate size distributions with a sharp transition region, one for the jammed case, and the other for the fluidized region. If the simulation would have been started from an initially fluidized state, the aggregate size distribution would be almost homogeneous over the gap.
Fig. 12 The aggregate size distribution (i is the number of monomers in an aggregate) at different times during flow start-up at an applied shear rate of a = 24 s−1 (see Fig. 8) for the superimposed aggregation kernel case. (a) Initial jammed configuration (b) distribution broken down close to the rotor (c) transient shear band has evolved (d) reached steady state with (trapped) shear band. |
A quite recent paper by Martin and Hu54 studies transient and steady-state shear banding for a thixotropic yield stress fluid (Laponite) in the shear rate controlled case. They observe three different transient shear banding scenarios depending on the amount of aging and the applied shear rate. At relatively low shear rates the final state is a shear banded one with the lowest shear rate of the sheared band being the critical shear rate. At relatively high shear rates the steady state was homogeneous. However, with sufficiently long aging times and intermediate shear rates the transient shear band was observed to be trapped. The behavior of the model here is quite similar and all of the different scenarios (except for an expanding low shear band) in particular trapped shear bands can be observed. The position of the trapped band is dependent on the initial state, which is here chosen by the steady states at very low shear rates while in54 the initial state is a function of aging time. However, the cases are equivalent in the sense that both longer aging time and lower shear rate imply a more jammed initial state. The trapped shear band scenario implies further that the simple lever rule,48 which states that the lowest shear rate of the flowing band is the intrinsic critical shear rate, does not always apply. However the model here does not predict any change in the critical shear rate. Contradicting experimental findings on this matter appear in the literature.54
In the model, determined by the active aggregate growth mechanism two distinctively different cases appear – the so-called simple and thixotropic fluids, characterized by monotonic and non-monotonic flow curves, respectively. Here, the non-monotonic flow curve of thixotropic fluids rises from a shear independent aggregation term originating from the Brownian movement of the particles, which dominates the aggregate growth at small shear rates, while its influence is negligible at high shear rates. As a result, at high shear rates, both fluid types behave the same, while at low shear rates additional complications appear in the thixotropic fluids. At zero shear rates, the simple fluids do not evolve in time, while the thixotropic fluids tend to age, following precisely their experimentally observed character.
Under non-homogeneous shear, the model fluids were shown to also follow qualitatively the recent experimental findings of spatially resolved rheometry. Start-up transients of the model fluids showed transient flow heterogeneities in good agreement with these experiments. The fluidization of both fluid types followed a similar power-law behavior at high shear rates nicely producing the experimental findings. Furthermore, steady state spatial heterogeneities, i.e. shear bands in thixotropic fluids and shear localization in simple fluids, suggested according to experiments to fundamentally categorize the two types, were reproduced by the present model, based on two distinct physical aggregate growth mechanisms. This suggests that from the aggregating unstable colloidal suspensions one can find both simple and thixotropic suspensions by modifying the particle size. Smaller particles have usually better diffusivity and, based on the present analysis, make the corresponding suspensions more thixotropic.
Based on experimentally known facts, the perikinetic aggregation is only relevant for aggregates below a few micrometers size. In general, larger particles have smaller diffusivities. In the aggregation kernel applied here, this feature is only taken into account by the prefactor, and the collision frequency is based on the aggregate size differences. Furthermore, the kernels are strictly valid only at the dilute limit, as they do not take into account the collective motion of the particles, and multiple particle collisions. Such improvements to the present kernels would be of considerable interest, but have turned out to be highly non-trivial. However, despite the fact that such modifications would constitute an important factor in the quantitative fitting of such models against experimental data, they can be expected to play only a minor role in the qualitative model behavior discussed in the present paper.
As demonstrated earlier by us50 the elastic stress loading can have important influence in the start-up behavior of complex fluids. Here, the elastic deformation of the colloid structures is not taken into account. Instead, we took the approach of a simple viscous constitutive relationship. Such a mechanism could accommodate for an elastic energy storage, which could be especially important during the start-up transients. Since in most applications complex fluids are in a transient state, finding a proper way to connect the elastic deformation of the aggregates to the fluids continuum elastic modulus will be an important future challenge.
This journal is © The Royal Society of Chemistry 2014 |