Rheology Dynamics of Aggregating Colloidal Suspensions

a 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.


Introduction
The rheology of complex uids shows both uidand solid-like characters due to their composition. Generally, these are mixtures of two phases with either gas, solid, or liquid of different viscosity, embedded in liquid. This embedded material can form structures depending on the external agitation and temperature. Under certain conditions, these structures can even ll the whole volume in an arrested conguration. 1 Then, macroscopically the uid appears to have load carrying capabilities, or in other words, a yield stress, associated to the solid-like state. 2 Furthermore, in such systems, when the external agitation is increased, a uid-like state is recovered. Even in the uid state, these materials usually show highly nonlinear behavior.
Unstable colloidal dispersions are one of the great examples of uids exhibiting many of the described phenomena. 3 These materials are vast in technological applications. Oen aggregation is intentionally induced due to the requirements of the process. For instance, water removal through ltration 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 lters. 4 Furthermore, since aggregation processes are usually reversible, 5,6 the colloids can be re-dispersed aer the required processing.
The complex rheology related to the change of the internal state in the suspensions is still rather poorly understood. Oen, the complex uid 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 uidity, 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 uids rheological properties. [8][9][10][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 ow can experience both, peri-and orthokinetic aggregation. Of these, one or the other dominates, depending on the ow intensity and particle properties. To reect 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 uid rheology are the Shear Transformation Zone (STZ) 21 and the So 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 so glassy materials formed of repulsively interacting so 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 uid is modeled via mesoscopic uid 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 uids 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 uids. We implement the model in a realistic measuring geometry, and compute the uid'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 uidization process.

Model
From rst principles, the aggregation dynamics is modeled in the form of mass balance equations originally written by Smoluchowski. 25 The discretized version of the mass balance equation reads k ðaÞ ði À j; jÞn iÀj n j À X N j¼1 k ðaÞ ði; jÞn i n j À k ðbÞ ðiÞn i þ X N j¼iþ1 bði; jÞk ðbÞ ð jÞn j ; where the subscripts i and j refer to different size classes, n i is the number of aggregates having i monomers per unit volume, and b(i, j) is the fragment size distribution 26 (here binary fragmentation). The kernels k (a) and k (b) depend on the type of ow and particle interactions. Only very selected cases have analytical kernels to begin with. Here, we apply the analytical aggregation kernels for simple laminar shear (orthokinetic aggregation), and thermal motion (perikinetic aggregation) to model the coagulation kinetics of Brownian particles under laminar shear ow. Superimposing the two processes gives the aggregation kernel where _ g is the effective shear rate, r is the aggregate radius, c o (c p ) is the collision efficiency for the orthokinetic (perikinetic) aggregation, T is the temperature (here T is xed to 300 K), and h 0 is the viscosity of the suspending uid. The collision efficiency terms correct the aggregation kernel, which are derived for neutral impermeable particles, and thus in practise overestimate the collision frequency if the particles exhibit nonneutral interaction potentials or have reduced permeabilities. 27 For the following analysis c o is assigned a constant value of 0.03 to account for this, whereas c p is varied to simulate the sensitivity of the aggregates to Brownian motion. In particular c p is assigned values from 0 to c o . If the model would be quantitatively tted against experimental data, these terms could be computed using either a capture cross-section method 27 for the orthokinetic kernel, or a stability ratio 14 for the perikinetic kernel and the terms would depend on the size and fractal dimension of the aggregates, and interaction potentials. Due to the lack of better alternatives, we use the kernels of eqn (2) irrespective of the concentration range, despite the fact that they are strictly speaking derived in the limit of low concentrations and ignore the collective motion of the particles and multiple particle collisions. All these effects can be expected to scale with the shear rate in the orthokinetic case and with the temperature in the perikinetic case. These scalings are similar to the original kernels, and thus these effects would change only the kernel prefactor, and would be irrelevant to this qualitative analysis. The perikinetic aggregation term introduces a negative slope to the ow curve at shear rates below the associated critical shear rate, _ g c . In particular the viscosity starts to diverge rapidly as the shear rate decreases below _ g c . A cut-off is applied to the perikinetic aggregation term at sufficiently high viscosities. This is done for two reasons, the rst one is due to obvious numerical issues, while the other one has a physical interpretation. Applying a cut-off that prevents the viscosity from diverging to innity is common in numerical simulations. 28,29 The applied cut-off can be interpreted as modeling the occurrence of trapping of the aggregates at high viscosities, which tends to slow down the aggregation rate. Due to this effect a suspension le at rest will age due to perikinetic aggregation approaching towards a steady state conguration as the rate of aggregation goes to zero due to jamming. Due to the lack of sophisticated models a simple tanh-type of cut-off was used throughout, where c p was scaled above a set viscosity Here h c is the cutoff viscosity for the perikinetic aggregation; here h c ¼ 40 Pas was used. The type of cutoff and the value of h c clearly have an effect on the resulting rheological behavior. However, based on extensive tests, the results below are not qualitatively affected by the cut-off function. 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 uid ow. 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 denite rst 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 aggregation 18 where F c is the aggregate effective bond force resisting the shear splitting. With this kernel the fragmentation rate increases rapidly when the hydrodynamic force (here estimated by Stokes drag) acting on an aggregate grows larger than the bond force. 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 reads 18,34 where h 0 is the uid viscosity, and k is a critical exponent indicating the divergence of the viscosity near the maximum packing fraction f max . The jammed volume fraction f is given by a summation over the partial volume fractions of the entire size distribution where the aggregates are assumed to cover the volume inside their rotational sphere, the radius of which is obtained using their fractal dimension (d F ) and the number of monomers (i) (with c, the lacunarity, set to unity here).
Neglecting the aggregate elastic stretching, assumes a viscous stress response s ¼ h _ g. For the sake of a realistic starting point we use parameters that were obtained earlier 18 by tting the orthokinetic model against rheological data for latex suspensions. 16 In particular 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 xed 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 f close to f max (see eqn (5)). We tested different number and spacing of size classes. These cause small systematic shis in the stress values. However, other quantities such as the local ow proles in the Couette are not inuenced 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 r 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 ow and no-slip boundary conditions, results in a function for the local shear rate 37 where U a (U b ) are the inner(outer) radial velocities, R b (R a ) the outer(inner) radii of the Couette and h(r) the local viscosity as given by eqn (5) with the corresponding f. Here we sample the radial shear rate distribution at least at 50 discrete positions. A separate homogeneous shear model is evolved simultaneously at each sampling point. During each iteration, the Stokes equation is updated to match the new structure. In what follows, this model will be called the heterogeneous shear rate model. When a ow curve of a uid 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 where R b (R a ) and U b (U a ) are the outer(inner) cylinder radii and angular velocity respectively. The simulations were carried out for a Couette of dimensions R a ¼ 2.39 cm and R b ¼ 2.50 cm, which are typical values for a narrow gap Couette and similar to the geometry used in the rst observations of transient shear banding in Carbopol gels. 40,41 In order to simplify the following discussion we dene the radial position as where d˛[0,1].

Homogeneous model
Computing the steady state stresses as a function of the applied shear rate with different kernels at various concentrations using the homogeneous shear model gives the ow curves plotted in Fig. 1. The concentration is here given in terms of which is the volume fraction with a pure monomer distribution, directly proportional to the solid mass fraction when the density is known. The ow curves show the typical shapes of shear thinning at small concentration, simple yield stress at increased concentration, and thixotropic uids when including the perikinetic aggregation term (c p > 0). Thixotropy is here understood as a property of the material to show aging. In order to simplify the discussion, we will later on refer c p ¼ 0 as the orthokinetic and c p > 0 as the superimposed case.
In the steady state ow curves the stresses, including the yield stress, drop with decreasing concentration, following experimental ndings of the behavior of unstable colloids. 16 In the model, the critical shear rate _ g c is a power-law function of both the concentration and the weight of the perikinetic aggregation, c p . The concentration dependence is observed also experimentally. 42 The dependence on perikinetic aggregation is experimentally realized by changing the particle interaction potential. 43 3.2 Heterogeneous model: steady state Now, we turn the discussion to the heterogeneous shear model. As mentioned before, classical rheometry assumes a constant shear rate throughout the gap, converting the angular velocity difference to apparent shear rate using relationships similar to eqn (8). On the other hand, eqn (7) shows that the local viscosity and shear rate are coupled. Owing to the fact that in a complex uid the viscosity is a function of the shear rate, it varies throughout the gap. Thus, the constant shear rate approximation has been considered useful only for sufficiently small gaps, 38 where the shear stress is almost constant. However, even in that case the use of eqn (8) (and any similar relationship that assumes a constant viscosity) can be questionable if the underlying uid has a highly non-linear ow curve, for instance close to the yield stress regime. Demonstrating such problems we plot in Fig. 2 the steady state ow-curves one would obtain when measuring a colloidal suspension in a Couette device using eqn (8) alongside with its real intrinsic behaviour as obtained using the homogeneous shear model. The ow curves shown in Fig. 2 are computed starting each sampling point from two different initial states. The initial aggregate size distributions are obtained from the homogeneous shear model either at _ g 0 ¼ 100 s À1 , corresponding to a uidized initial state, or at _ g 0 ¼ 10 À5 s À1 , corresponding to a jammed state (further on referred as uidized and jammed initial states, respectively)thus this is not a sweep simulation per se, but the steady states  starting from the different initial distributions correspond to the ones that would be obtained by doing down sweep and up sweep shear ramps, respectively. The deviations from the intrinsic behaviour are easily spotted. Firstly, the measured ow curves have somewhat higher stresses than the intrinsic ones. Secondly, no negative slope is reported by the device. In fact from these shear rate controlled apparent steady state results it is not possible to see whether a uid shows yield stress behaviour or a diverging stress at some critical shear rate. Stress controlled runs can result in better results in this regard. 44 Moreover, in the superimposed case the reported ow curve is dependent on the initial state, whereas for the case with a monotonous ow curve, the one reported using the Couette device is unique. Non-unique ows have also been experimentally observed in thixotropic systems (e.g. ref. 3 and 45).
The cause of the deviations from the intrinsic ow curve is the heterogeneous shear prole, 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 uid, 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 proles. When the uid is initialized in a state having smaller viscosity than the steady state one, the static shear band has the width (h SB ) determined by _ g c in accordance to the simple lever rule 48 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 uid'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 ow curves. Steady state radial shear rate proles 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 rst scenario is the case where the uid is initially in a uidized 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 ow curve. The second scenario is shear rate controlled and starts from a jammed state shown in Fig. 4(b). As mentioned before, in that  Local stress-shear rate values at steady state with the superimposed aggregation kernel at different applied engineering shear rates ( _ g a ¼ 0.05, 4, 24, 100 s À1 ). When initial state is (a) fluidized _ g 0 ¼ 100 s À1 (b) jammed _ g 0 ¼ 10 À5 s À1 . case the resulting shear rate prole 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 uidized state gives homogeneous shear proles 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 conguration. For such situations the local shear rate-stress pairs are plotted in Fig. 5(b). Now, the le, low stress branch of the ow 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.

Heterogeneous model: transient behaviour
The response of the orthokinetic case to shear transients resembles the behavior of simple structural models. 46,50 Starting from an initial state corresponding to high shear rates, while imposing a relatively moderate shear rate, shows the growth of shear localization as plotted in Fig. 6. There, as a response to the enhanced growth of the structure at the outer cylinder due to the lower shear rate, the shear rate decreases further, until in the steady state, part of the gap exhibits a very small shear rate. The transition from this quiescent region to the highly sheared region is in this case continuous. The nal conguration shown in the plot is the steady state corresponding to _ g a ¼ 0.7 s À1 . The stress relaxation, shown in the upper panel, occurs in an expected manner, gradually increasing towards the steady state.
Shear rate transients from jammed states to uidized ones are shown for the orthokinetic case in Fig. 7 and for the superimposed kernel in 8. In both cases, the ow rst localizes at the rotor aer 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 uidization occurs through a transient shear banding phase, nally ending in a linear velocity prole 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 rst 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 uidized 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 uids where the structural changes between the ow and no ow states are signicant. Divoux et al. addressed the issue of uidization experimentally for the case of a simple yield stress uid (i.e. not thixotropic) using extensive rheometry combined with local velocimetry measurements. 40,41 The main conclusions drawn are that the uidization time (reaching a homogeneous ow prole) can take a considerably long time. Moreover, they found that the uidization time follows power law scalings, with different exponents for the shear rate and stress controlled modes. In the present model, the uidization time, here dened 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 uidization time for shear rate controlled runs as a function of the engineering shear rate is shown. In the orthokinetic case, the uidization time can be well tted with a power law having the exponent a z 1. Changing the initial conguration (to more or less jammed) affects only the prefactor but not the power-law exponent.
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 conguration has a static shear band, the uidization time is not well dened (of course a timescale to reach a steady state can be dened, however). For intermediate shear rates, for which the steady state is not shear banded, the uidization time follows a powerlaw with an exponent that increases with c p (see Fig. 10). At very high shear rates ( _ g [ 100 s À1 ) the orthokinetic and superimposed cases coincide.
Next the focus is turned to the uidization process with stress controlled runs. The uidization times as a function of the reduced stress are shown in Fig. 11. Here, following the denitions used in the experiments, 40 the reduced stress means simply subtracting the yield stress from the imposed stress values. This way, the uidization times can be tted by a power-   law with the exponent b z 1.6. However, the exponent is quite sensitive to the used yield stress value (here s c ¼ 1.519 Pa obtained from the intrinsic ow curve has been used). At very high stresses the uidization time scales with an exponent b z 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 uidization time can be dened. Changing the initial conguration (to more or less jammed) does merely affect the pre-factor but not the exponent. Also the scaling at high stresses with an exponent b z 1 and the scaling at the low stress limit with b z 1.6 are independent of the concentration as in the shear rate controlled simulations.
To conclude the uidization times for the shear rate and stress controlled cases follow power law like behaviour. The yield stress can be obtained from the stress controlled uidization times. Divoux et al. observed that the ratio of the two uidization exponents gives the exponent of the Herschel-Bulkley model (s ¼ s o + h _ g a ) t to the measured materials (apparent) ow curve. 41 In the current model the value of the Herschel-Bulkley exponent scales with the k-exponent in the viscosity function eqn (5). However, the uidization time exponents were insensitive to changes in k. 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 gure 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 ow 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 uidized region. If the simulation would have been started from an initially uidized state, the aggregate size distribution would be almost homogeneous over the gap.

Comparison to experiments
With the model utilized here several aspects of yield stress uids have been investigated. In particular the change from simple yield stress to thixotropic yield stress uids and the different rheological phenomena like shear banding associated with these have been addressed. Next we compare available experimental ndings with the model at hand. Changing the amount of perikinetic aggregation changes the simple yield stress uid to a thixotropic one. Tuning a system thixotropic has been experimentally addressed several times. 43,52,53 In particular Ragouilliaux et al. 43 observed that at high shear rates the ow curves for both cases, simple and thixotropic, are similar. This agrees qualitatively with what is observed in the model here. The origin of the thixotropy has also been addressed several times. In ref. 42 Coussot et al. show that by considering the ratio of the characteristic relaxation time and the restructuring time the steady state shear banding can be understood. In the model here the origin of the rheological features can be directly traced back to changes in the microstructure. The shear rejuvenation and aging observed in thixotropic yield stress uids 12 can here be understood by the shi of the aggregate size distribution towards smaller aggregates due to fragmentation and due to aggregate growth resulting from orthokinetic and perikinetic aggregation respectively.
A quite recent paper by Martin and Hu 54 studies transient and steady-state shear banding for a thixotropic yield stress uid (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 nal 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 in 54 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 owing 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 ndings on this matter appear in the literature. 54

Conclusions
We studied a continuum colloidal model showing ow behavior of typical complex uids. The rheological model based on the population balance equation is two-way coupled to the Stokes equation of a Couette geometry. We found this model to reproduce and explain many of the features familiar from these experiments in both simple yield stress uids, as well as in the thixotropic ones. The results obtained here further reinforce conclusions about static and transient shear banding of earlier modeling work. 46,50,51,55 Further the transition from simple to thixotropic yield stress uids is here a continuous one, in agreement with recent papers. 42,49 However due to the model being purely viscous, effects related to elasticity such as stress overshoots, are not present. In addition the transient behaviour of thixotropic yield stress uids, has been explicitly addressed from the modeling perspective. The corresponding results are well in line with experimental results. 54 With the model, the coupling between the microstructure (aggregate size distribution) and the rheological properties could easily be identied, information that is not obtained from simple structural models.
In the model, determined by the active aggregate growth mechanism two distinctively different cases appearthe socalled simple and thixotropic uids, characterized by monotonic and non-monotonic ow curves, respectively. Here, the non-monotonic ow curve of thixotropic uids 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 inuence is negligible at high shear rates. As a result, at high shear rates, both uid types behave the same, while at low shear rates additional complications appear in the thixotropic uids. At zero shear rates, the simple uids do not evolve in time, while the thixotropic uids tend to age, following precisely their experimentally observed character.
Under non-homogeneous shear, the model uids were shown to also follow qualitatively the recent experimental ndings of spatially resolved rheometry. Start-up transients of the model uids showed transient ow heterogeneities in good agreement with these experiments. The uidization of both uid types followed a similar power-law behavior at high shear rates nicely producing the experimental ndings. Furthermore, steady state spatial heterogeneities, i.e. shear bands in thixotropic uids and shear localization in simple uids, 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 nd 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 modications would constitute an important factor in the quantitative tting 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 us 50 the elastic stress loading can have important inuence in the start-up behavior of complex uids. 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 uids are in a transient state, nding a proper way to connect the elastic deformation of the aggregates to the uids continuum elastic modulus will be an important future challenge.