3D simulations of the flow of thixotropic fluids, in large-gap Couette and vane-cup geometries

J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Contents lists available at ScienceDirect

Journal of Non-Newtonian Fluid Mechanics journal homepage: www.elsevier.com/locate/jnnfm

3D simulations of the flow of thixotropic fluids, in large-gap Couette and vane-cup geometries Andrei Potanin Colgate Palmolive, 909 River Rd., Piscataway, NJ 08855, United States

a r t i c l e

i n f o

Article history: Received 26 February 2009 Received in revised form 13 October 2009 Accepted 7 January 2010

Keywords: Vane Thixotropy CFD

a b s t r a c t Models of the vane-cup and Couette rheometers are compared using computational fluid dynamics as well as approximate solutions. Thixotropy of the fluid is incorporated by means of a model based on experimental data for various toothpastes. Parameters of the model are calculated by fitting the results of step-shear tests in Couette geometry and are subsequently used to predict torque for similar tests in vane-cup geometry. Calculations for the thixotropic model were compared to the calculations for a timeindependent equilibrium model with the same steady-shear rheology. There are two important findings in this work. First, at constant angular velocity, torque in the vane-cup rheometer for the thixotropic model turns out greater than for the equilibrium model, because the structure of the thixotropic fluid has not had enough time to reach equilibrium. This has implications for rheometry and the modeling of mixing operations. In both cases, torques may be underestimated if the standard equilibrium model is used in calculations. The second finding is relevant to rheometry. As is well known, equilibrium model predicts a lower torque value for a vane-cup geometry than for an equivalent Couette geometry. We found that taking thixotropy into account either makes the difference less pronounced or in some cases actually makes torque for vane-cup higher than for Couette. End effects are also analyzed. © 2010 Elsevier B.V. All rights reserved.

1. Introduction Despite a great number of modeling and experimental investigations of thixotropy [1,2], there is a certain gap between a practical industrial approach and fundamental rheological studies. Rheologists prefer to deal with simple geometries with well-defined shear rate, such as cone-plate, whereas practical calculations require analysis of complicated geometries. Commercial computational fluid dynamics (CFD) packages, such as Fluent or Flow-3D, which are becoming widely used in industry for such calculations, provide only limited or no standard options for thixotropic effects. Most of the time thixotropic effects are completely neglected in CFD. For instance, mixing operations of pastes are usually analyzed in terms of simple yield-stress models, such as Herschel–Bulkley or Casson. Since practically all yield-stress materials are actually thixotropic, the local viscosity depends on the extent to which structure has reached equilibrium. Very often the structure has no time to be broken down between blades, which means that neglecting thixotropy leads to an underestimation of the torque on the shaft of a mixer. A related problem arises in rheological measurements. Many testing protocols for concentrated suspensions and pastes use nonconventional rotational geometries with no simple mathematical

E-mail address: andrei [email protected]. 0377-0257/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.jnnfm.2010.01.004

relation between the angular velocity and local shear rates in the sample. The problem is even more serious in transient shear experiments, since shear rate is not only heterogeneous, but also time-dependent in ways which are impossible to predict without some sort of a rheological model. Vane-cup geometry provides a relevant example, since it is now becoming a standard tool in many industrial applications [3], although the data analysis is often inadequate, especially for time-dependent flows. Mixing and vane-cup rheometry are closely related. Vane-cup itself may be thought of as the simplest model of a mixer. An approximation is often made in which vane is considered as a solid cylinder encompassing its blades [4]. Some 2D computer simulations were performed to demonstrate that strongly shear-thinning fluids tend to circulate in the zone between blades without leaving this zone. Thus, the flow around a vane is indeed similar to the flow around its encompassing cylinder [5]. More recently, 3D calculations have basically confirmed this [6]. Thixotropy may affect these assumptions in the same way it affects flow in any mixer since the structure may not have enough time to reach equilibrium between blades. Large-gap Couette flow between two infinite cylinders in transient regimes was modeled by Roussel et al. [7]. They calculated velocity distribution of a thixotropic fluid in this geometry by extending a well-known approach of integrating shear rate over shear stress in the gap. Such calculations are possible only after a rheological model is selected. Hence, in order to be usable, the

300

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Nomenclature f () ˙

function which defines explicit shear rate dependence of thixotropic viscosity in Eq. (5) L, Ls length of internal cylinder (or vane) and shaft k parameter of f (), ˙ Eq. (6) K number of discretization points in Eq. (18) mr , mb parameters of tc , Eq. (9) n parameter in Casson equation nc , nr , nb exponent in Eq. (8) and its values for recovery and breakup in Eq. (9) p pressure R, R0 , Rs radii of internal cylinder (or vane), cup and shaft r radial coordinate T torque on the shaft of the cylinder (or vane) Tr , Tb recovery and breakup times at reference shear rate, Eq. (9) t time tc characteristic time of recovery or breakup, Eq. (8) m m ˙ r, tbreak ≡ Tb (˙ ref /) ˙ b recovery and trec ≡ Tr (˙ ref /) breakup times, Eq. (9) trot ≡ /2˝ the time taken for the vane blades to rotate 90◦ v fluid velocity x ≡ r/R dimensionless radial coordinate discretization points in Eq. (17) xj , ˙ ˙ rate of strain tensor and shear rate ˙ 1 , ˙ 2 initial and final shear rates in a step-shear test parameter of f () ˙ ˙ 0 2 ˙ app ≡ 2˝/(1 − (R/R0 ) ) apparent shear rate on the surface of an infinite cylinder in a Newtonian fluid ˙ ref reference shear rate, Eq. (9) ˙ min cut-off shear rate in Eq. (10) ␦ unit tensor ε ε ≡ Rs /R  relaxation time , eq , ∞ viscosity, equilibrium viscosity and high-shear viscosity (parameter in Casson equation) ,  eq structural parameter and its equilibrium value  in initial value of the structural parameters in CFD  density ␴, , y stress tensor, shear stress and yield stress (parameter in Casson equation) initial and final shear stress in a step-shear test 1, 2 min minimum shear stress in an undershoot wall shear stress at the wall of a cylinder ˝ angular velocity ˝1 , ˝2 initial and final angular velocity values in a stepshear test

model has to be flexible to include various types of thixotropic behavior encountered in pastes. However, such flexibility requires use of a large number of fitting parameters. This paper is written by an industrial rheologist who has experienced first hand a gap between rheological and CFD communities with regards to modeling of thixotropic materials and intends to demonstrate how even a relatively simple thixotropic model can be used to bridge this gap. The goal is not to formulate another thixotropic model, but to show how to improve modeling accuracy of mixing of thixotropic materials in CFD. From this standpoint, Couette and vane-cup geometries may be thought of as but two simple mixers. How effective they are as mixers is irrelevant here, since we are only concerned with the torque on their shafts. Of these two geometries Couette has an advantage in that it allows an iterative solution in integrals (although some CFD modeling is

needed to account for end effects) whereas vane-cup can only be analyzed by CFD. Thus, this paper is organized as follows: we first analyze Couette data to determine model parameters and then use these parameters to calculate flow in vane-cup. The way one reads this paper will depend on one’s background and point of view: one may either consider vane-cup as just a mixer on which novel CFD approach is tested or as a rheometrical tool which is analyzed using CFD technique. 2. Experiments The results presented in this work were obtained for three toothpastes which we denote as Pastes 1, 2 and 3. These pastes were chosen simply because they represented different types of thixotropic behavior. The main differences between the pastes are in the type of abrasive and polymers they used. Pastes 1 and 2 contained precipitated calcium carbonate abrasive particles and Paste 3 contained silica particles. Paste 1 contained cellulose and xanthan gum, while Pastes 2 and 3 contained cellulose gum and carrageenan. Paste 3 also differed from the other two in that it contained a high molecular weight polymer component. A Malvern Bohlin HR-nano rheometer was used for all tests. Two geometries were used, large-gap Couette with a serrated internal cylinder and 4-bladed vane in the same cup. Both were made from stainless steel. A sketch of the vane-cup geometry is shown in Fig. 1. The cylinder for Couette geometry had 35 vertical serrations 0.5 mm deep to prevent slippage. These geometries are similar to those supplied by Malvern except that they had thinner shafts (3 mm in diameter) so as to reduce their contribution to torques and the cylinder had flat bottom so as to exactly encompass the vane. The Malvern rheometer internal cylinder (or vane) rotates with angular velocity ˝. For convenience in what follows, we will use the so-called “apparent shear rate” defined as the shear rate on the surface of an infinite cylinder in a Newtonian fluid, ˙ app = 2 /(1 −

Fig. 1. Schematic representation of the vane-cup geometry. Couette geometry had the same dimensions with the internal cylinder encompassing the vane. All dimensions are in millimeters.

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312 2

(R/R0 ) ), where R and R0 are radii of the internal cylinder (or vane) and the cup. Accordingly, apparent shear stress will be defined as the shear stress on the wall of that imaginary cylinder, T/2R2 L, where T is the torque on the shaft of either the cylinder or the vane and L is the length of the vane (cylinder). All experiments reported in this work were transient step-shear tests. Each test was performed on fresh samples according to the following protocol. First a pre-shear at the rate of 30 s−1 was applied for 1000 s to reach an equilibrium reproducible torque reading. The recovery tests were performed by switching the shear rate from its pre-shear value to the final value, which ranged from 0.03 to 1 s−1 . The breakup tests, on the other hand, were always preceded by a breakup test from 30 to 0.1 s−1 . Then, after shearing at 0.1 s−1 for 3000 s, the shear rate was switched to its final value between 1 and 30 s−1 . This step-shear protocol is widely used today for characterization of thixotropic materials, even though it is sometimes criticized as too simple (or even old-fashioned) by those who favor more dynamic protocols with continuous increase or decrease of angular velocity. We decided to choose this step-shear protocol mainly because in practical industrial rheology this approach still dominates and probably will dominate in the nearest future. We were also motivated by a recent series of papers by Dulaert and Mewis [2] who used this same protocol on a model system. 3. Theory In this work, we use various models to solve different problems in two different geometries, Couette (or concentric cylinders) and vane-cup. We apply CFD to study both of them. For vane-cup this requires 3D modeling, while for Couette the problem is reduced to 2D due to axial symmetry. Also, we can get a simplified solution of the Couette problem by reducing it to 1D problem of infinite cylinders. For each of these geometries there are two different problems to consider, one dealing with variable and another with constant angular velocity. In fact, this work deals with only one variablevelocity problem which is the step-velocity (or step-shear) problem with angular velocity (or apparent shear rate) changed by step as described in the previous section. In what follows we will consider two models, equilibrium and transient. An equilibrium model postulates a certain timeindependent viscosity being a function of local shear rate, whereas a transient model introduces time dependence of rheological properties. In this work we consider a simple thixotropic model, in which at any constant shear rate in the long-time limit a unique equilibrium solution exists. For this model in Couette geometry constant-velocity problems may be thought of as a particular case of the variable-velocity problem with the given angular velocity maintained for an infinitely long time. In vane-cup geometries, even in constant-velocity problems and even for our simple thixotropic model, local shear rates vary, so that equilibrium and transient models give different solutions. We will compare these solutions in the next section. But let us first define the models.

301

where , v, ␴, ␦, and p are fluid density, velocity, stress tensor, unit tensor and pressure, respectively. Viscosity is a function of shear rate, , ˙ and structural parameter, : ˙  = (, ) ˙ .

(3)

To complete the model we need to specify the viscosity func˙ if a tion, (, ). ˙ It should attain the equilibrium viscosity, eq (), constant shear rate is maintained. Casson equation provides convenient flexibility and fits equilibrium data for various pastes:

 eq () ˙ ≡ n∞ +

 n 1/n y

˙

,

(4)

where ∞ , y and n are high-shear viscosity, yield stress and an exponent which defines the shape of the equilibrium flow curve. Equilibrium model assumes that viscosity is always set to its equilibrium value. To define our thixotropic model, we need to specify the dependence of viscosity on the structural parameter, . There are two commonly used approaches. The first approach, which may be traced back at least to Moore [8] and is still widely used [7,9], is to assume that  depends only on , i.e., there is no explicit dependence of viscosity on the current shear rate. Another approach (suggested in a different context by Tiu and Boger [10] and De Kee et al. [11]) is to assume that viscosity depends linearly on both the structural parameter and equilibrium viscosity, i.e.,  = eq (). ˙ Some more complex models [12] reduce to this latter approach in the low-shear limit, which is relevant here. It turns out that neither approach quite complies with our step-shear data. To explain this predicament, consider a down-shear step in which the shear rate changes from a high-shear rate, ˙ 1 , to a lowshear rate, ˙ 2 . Correspondingly, structural parameter changes from its initial equilibrium value, eq (˙ 1 ), to a non-equilibrium value, , which eventually approaches a new equilibrium, eq (˙ 2 ). Shear stress goes from its initial value, 1 , through an undershoot with its minimum, min , approaching a new equilibrium, 2 . Its undershoot value depends upon the initial value of the structural parameter, eq (˙ 1 ), since the structure has not recovered yet. Had viscosity been dependent only on  (the first approach), we would have had min ∼eq (˙ 1 )˙ 2 . On the other hand, had viscosity been proportional to its equilibrium value (the second approach), we would have had min ∼eq (˙ 1 )eq (˙ 2 )˙ 2 . Since in the low-shear limit eq () ˙ ˙ in the first case min would be linear in ˙ 2 , while approaches y /, in the second case it would not depend on ˙ 2 at all. Experiments suggest something in between, namely that min is proportional to ˙ 21−k , where k usually ranges between 0.4 and 0.7. Thus, still retaining a linear dependence on the structural parameter, we generalize the formula for viscosity as follows: (, ) ˙ = f (), ˙

(5)

where f () ˙ is introduced to account for the dependence on the current shear rate for which the simplest formula is



f () ˙ = ∞ 1 +

 ˙ k  0

˙

,

(6)

where ˙ 0 and k are parameters, and ˙ eq () . f () ˙

3.1. Models

˙ = eq ()

The usual mass and momentum conservation equations for isotropic, incompressible and isothermal fluid in the frame of reference rotating at constant angular velocity ˝ are as follows:

To describe kinetic effects we need to define the transport equation for the structural parameter. Following what is sometimes referred to as an “indirect microstructural” approach [1], we postulate the following nc th order kinetic equation for the structural parameter:

∇ · v = 0, ∂v  + v · ∇ v + 2 × v +  × r × r = ∇ · (−pı + ␴), ∂t

(1) (2)



(7)

nc

 − eq () ˙  ∂ + v · ∇ = ± tc () ˙ ∂t

,

(8)

302

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

where tc () ˙ is the characteristic time of the recovery or breakup processes to both of which Eq. (8) applies with plus or minus sign, respectively:



tc = trec ≡ Tr

˙ ref

˙ ref

for recovery,

( < eq )



, nc = nb

 −˙ 

˙ min

for breakup.

( > eq ) (9.2)

,

(10)

which effectively cuts viscosity above y /˙ min . In what follows, except when otherwise stated, calculations were performed with ˙ min = 10−4 s−1 . Although, microphysical interpretation of the constitutive equations introduced above is not a subject of this work, some remarks on Eqs. (5) and (6) may be made to put them into a perspective of the existing models. Indeed, these equations may be rationalized by separating shear stress into some “viscous” and “elastic” (or “aggregates” and “networks”) parts as many authors have done recently for dispersions (e.g., [15,16] or [17]). The “viscous” part, ∞ , ˙ may be thought to represent slowly growing aggregates and is therefore directly related to the structural parameter which presumably depends on their size. The remaining “elastic” part represents a network. 3.2. 1D solution for simple shear In simple shear Eqs. (3) and (8) reduce to = (, ) ˙ , ˙ d =± dt

 nc  − eq () ˙  tc () ˙

(11) ,

eq (˙ 1 ) eq (˙ 2 )



e−t/tc .

(13)

3.3. 1D solution for Couette flow (infinite concentric cylinders)

mb

˙

1−

(9.1)

where ˙ ref is a reference shear rate value, which we set to 1 s−1 . Similar power laws have been often suggested for the recovery and breakup functions [2]. Eq. (8) also is often used in thixotropic modeling. If the exponent nc is set to 1, it may be rewritten in somewhat more appealing terms of the first-order kinetics of bonds creation/breakup (“direct microstructural” approach). However, experiments described below suggest qualitatively different kinetics on recovery (in down-shear steps) and breakup (in up-shear steps) as well as different kinetics for different pastes. It is worthwhile to note that the above thixotropic model may be thought of as a generalization of the model commercially implemented by Flow Science Inc. in their Flow3D code [13]. Our model reduces to theirs, if we put f = ∞ in (5) and nc = 1 in (8), so that viscosity becomes proportional to the structural parameter and kinetic equations are simplified to the first order. While our approach increases the number of model parameters, this turns out to be necessary to quantitatively describe experimental data and to correctly model industrial processes. An important remark has to be made regarding the numeric implementation of this model. For the actual calculations viscosity needs to remain finite in the low-shear limit. This is usually achieved by cutting-off its divergence below certain shear rate, ˙ min , which is selected low enough not to affect relevant calculations. To this end, following Papanastasiou approach for Bingham plastics [14], we make the following substitution for the yield stress for all numerical calculations: y → y 1 − exp



= 2 − 2 , n c = nr

˙

 tc = tbreak ≡ Tb

mr

to low, ˙ 2 . Assuming for simplicity that kinetic equation for the structural parameter is of the first order, nc = 1, we get

(12)

where viscosity function is defined by Eq. (5). As a demonstration of how this model works, let us again consider a step test in which shear rate is decreased from high, ˙ 1 ,

In Couette flow, as long as infinite cylinders approximation applies, equations in cylindrical coordinates are similar to (11) and (12), except that the shear rate, ˙ = x∂˝/∂x, and other variables depend on the radial coordinate x ≡ r/R. Hence [18]:



R0 /R

˝= 1

(x) =

(x) ˙ dx, x

(14)

wall , x2

(15)

where wall ≡ (1) = T/2R2 L is the shear stress on the internal cylinder wall. For steady-shear the problem is further simplified since the function () ˙ is known, e.g., given by the Casson Eq. (4) as eq () ˙ , ˙ where fitting parameters are to be estimated by fitting T(˝) to the experimental data. For transient flow the solution is easy to obtain iteratively. We have wall (t) = x2 (x)( ˙ (x), ˙ (x, t))

(16)

At time t = 0, immediately after the step in shear, the structural parameter is known, (x, 0) = eq (˙ eq (x)), where ˙ eq (x) is calculated by solving equilibrium Eqs. (14) and (15) at ˝ = ˝1 . Hence, Eq. (16) gives (x) ˙ with the only unknown parameter, wall . Substituting (x) ˙ into (14) at ˝ = ˝2 we solve it to calculate wall (alternatively, one can differentiate (16) over x and solve the second-order equation for ˝(x) with boundary conditions ˝(1) = 0 and ˝(R0 /R) = ˝2 , as was done by Baravian et al. [19]). Subsequent time steps require only direct integration by Runge-Kutta method of the following equations: d(xj , t) dt

=±

 n (xj , t) − eq ((x ˙ j , t)) tc ((x ˙ j , t))

,

j = 1..K,

(17)

for K points into which the gap is discretized. Each time the local shear rate, (x ˙ j , t), is calculated from Eq. (16) and wall is recalculated by the following integral: wall =

R0 /R 1

∞ ˝

(18)

dx/x3 ((x, ˙ t), (x, t))

3.4. 3D solutions with Fluent Fluent fluid dynamics software commercial packages (versions 6.3.26 and 12.0.7 by Ansys) were used for calculations. Meshes were generated using the Gambit 2.4.6 package. Since vane is symmetrical with respect to 90◦ rotation, only one sector of the whole geometry needs to be simulated for vane-cup geometry and periodic boundary conditions are therefore used with opposite faces appropriately linked. Cross-section of the mesh is shown in Fig. 2. This mesh has 116,732 hexahedral cells and was generated by separating the whole fluid volume into sections and using the so-called Cooper scheme which projects a face mesh from one end of each volume to another. Our mesh for Couette geometry looked similar to vertical cross-sections in Fig. 2 but also used finer meshing (up to 5745 against 3008 in Fig. 2). In view of its cylindrical symmetry, the solution for the finite-length Couete problem was actually obtained in 2D-axisymmetric mode, which is provided in Fluent. As usual, we varied mesh size and shape in both cases so as to make sure that the solution does not depend on it.

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

303

but as we will show our results practically do not depend on it, because its effect is eliminated during the pre-shear. In other words, a pre-shear is as necessary for calculations as it is for the actual experiments and for the same reason, to condition the sample. 4. Results 4.1. Experiments

Fig. 2. 3D mesh.

As explained in the beginning of this section, we consider here equilibrium and thixotropic models. Our equilibrium model uses Eq. (4) and is implemented via a user-defined function for material properties written in C as customary for Fluent applications. Our thixotropic model was implemented via two user-defined functions. One function defined the evolution of the structural parameter, , according to the transport Eq. (8). Another function defined viscosity relating it to the shear rate and . The transport equation was solved together with flow equations using coupled algorithm and under-relaxation factors so as to achieve stable converging solution with tolerance at least 10−6 . Fluent solves the complete flow equations (2), although in this particular application inertial terms can be as well neglected since we are only dealing with low Reynolds numbers. Also, Fluent allows solving the problem in a moving reference frame, which corresponds to the actual experimental setup with rotating vane or internal cylinder or in a stationary frame with angular velocity defined on the cup. Again, at low Reynolds numbers both approaches give the same solution, although convergence is slightly better in the latter case. Torque can be estimated either on the cup or vane or internal cylinder. If the mesh is sufficiently fine, both estimates coincide, although generally the estimate on the cup tends to be more accurate. Two different Fluent solvers were used: steady and transient. The steady solver was used for the equilibrium model, and transient for the thixotropic model. As expected, the solution for the thixotropic model approaches the solution for the equilibrium model for Couette geometry in the long-time limit. However, in vane-cup geometry calculations with the thixotropic model significantly differ from those with the equilibrium model even in the long-time limit, as will be shown below. Transient calculations require initial conditions to be specified. Except where otherwise indicated, the structural parameter  was set initially at  in = 2. This may seem a rather arbitrary choice,

All experimental data reported here were obtained in stepshear tests performed as described in Section 2. Since each test is preceded by a long pre-shear, pre-conditions are thereby mostly eliminated. There are indeed stagnant areas in both geometries (between blades and on top and bottom) but, as we will show below, they do not contribute to the torque. Equilibrium flow curves are plotted in Fig. 3 for all three pastes. These data were collected in Couette geometry. Each point of each curve corresponds to a fresh sample which underwent a sequence of shear rate steps defined in Section 2. For Paste 1 these steps are shown in Figs. 4 and 5. Thus, shear stresses at shear rates of 1 s−1 and below were collected at the end of the recovery tests, like those in Fig. 4, each of which, in its turn, was preceded by a pre-shear at 30 s−1 . On the other hand, shear stresses at shear rates above 1 s−1 were collected at the end of breakup tests, like those in Fig. 5, each of which was preceded by recovery at 0.1 s−1 . The meaning of the theoretical curves on these figures will be explained in the next section. Note that we plotted the data both in logarithmic and linear coordinates so as to better illustrate our data and calculations at different time scales. Since, in what follows, long-time limits are especially important, all other data will be plotted only in logarithmic coordinates. Figs. 6–8 summarize all step-shear tests for the three pastes performed in Couette and vane-cup geometries. Plots (a) on the left correspond to recovery tests, while plots (b) on the right to breakup tests. Vane-cup data generally fall below Couette data, except at low-shear rates for Paste 1. Again, curves on these plots will be explained in the next section. 4.2. Calculations It is often assumed that at low-shear rates torque on the Couette cylinder can be estimated by integrating yield stress over its whole surface [3,4]. This may or may not include the shaft which is also a cylinder with radius Rs and immersed length Ls . In our case the shaft

Fig. 3. Equilibrium flow curves measured with Couette geometry. Curves represent fit with Casson equation.

304

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Fig. 4. Step-shear recovery tests with Couette geometry for Paste 1 with apparent shear rate changed as shown on the plot (in s−1 ). Solid lines, CFD calculations in axisymmetric geometry. Dashed lines, calculations for 1D geometry. For parameters see Table 1. Two plots are identical, except for scale, and shown so as to reveal both short- and long-time transient behavior.

Fig. 5. Same as in Fig. 4, but for breakup.

Fig. 6. Step-shear tests with both Couette (filled symbols, thick curves) and vane-cup (empty symbols, thin curves) geometries for Paste 1: (a) recovery and (b) breakup. Apparent shear rates indicated on the plots. Dashed lines on (b) correspond to the equilibrium model for vane-cup geometry.

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

305

Fig. 7. Same as Fig. 6, but for Paste 2.

was much thinner than the main cylinder (see Fig. 1): R = 7 mm, L = 30 mm, Rs = 1.5 mm, Ls = 24 mm. Simple calculation gives the ratio of the torque calculated by integration over the whole surface, including the shaft, to the torque calculated by integration over side surface only end-effect correction factor = 1 +

2R Ls + ε2 + O(ε3 ), 3L L

(19)

where ε ≡ Rs /R  1 (see [20] for estimates of the higher order terms). Only the first two terms are often retained. This gives a well know no-shaft correction also used for vane by Nguyen–Boger [4] and others. While the O(ε3 ) in (19) gives less than 1% contribution,

the O(ε2 ) term may be a source of some concern. However, the validity of this estimate becomes somewhat questionable as the shear rate increases, since even at very low-shear rates the flow around the shaft becomes mostly screened by the flow generated by the large cylinder so that the true contribution from the shaft is actually much smaller. In Appendix A we compare our 3D calculations to calculations for infinite cylinders and infinite vane-cup geometries to estimate end effects for various set of parameters. In the most relevant range between 0.01 and 10 s−1 , the error of Nguyen–Boger correction for Couette does not exceed 2%. The fact that this correction works quite well can as well be seen for Paste 1 in Figs. 4 and 5 which compare Fluent calculations to iterative calcu-

Fig. 8. Same as Figs. 6 and 7, but for Paste 3.

306

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Table 1 Fitting parameters for the three pastes.

1 2 3

y (Pa)

∞ (Pa s)

n

˙ 0 (s−1 )

k

Tr (s)

Tb (s)

mr

mb

nr

nb

74 58 24

3.2 1.2 4.7

0.7 0.4 0.3

70 30 30

0.6 0.7 0.5

40 5 20

15 35 40

0.6 1 1.6

1 1 1

1 1.5 3

3 3 3

lations in 1D described in Section 3.3 for the same set of parameters. Hence, in what follows all 1D Couette calculations will be corrected according to Nguyen-Boger correction, i.e., by multiplying torque or apparent shear stress by the factor 1 + 2R/3L. All calculations were performed with parameters defined for each paste in Table 1. These parameters were determined by comparing measurements in Couette geometry to the calculations performed according to the transient model defined in Section 3 in either 3D using Fluent or in 1D using iterative procedure of Section 3.3. Since the latter is much easier and sufficiently accurate, we used 1D geometry to do all the fitting and after that we did calculations in all other geometries with the same parameters. Having 12 parameters makes the model very flexible so that data for practically any paste can be fitted. In practice some of the parameters may be fixed. Thus, in what follows we considered breakup parameters, mb and nb , fixed at 1 and 3, respectively, which is good enough to describe available data. If one wants to think in microscopic terms, one may consider having mb = 1 as an indication that breakup is shear-controlled. High nb indicates that breakup curves are “stretched” or slow decaying as compared to recovery curves (except for Paste 3, in which all curves are equally “stretched”). As for other parameters, the following procedure was used. To start with we fitted equilibrium flow curves in Fig. 3 to determine the three parameters of the Casson equation (the first three columns in Table 1). Next, we fitted recovery curves, like those in Fig. 4. This allows for an estimate ˙ 0 , k, Tr , mr and nr . The first two parameters determine the depth of the undershoot on the curve, the next two define breakup time dependence on shear rate, while the last one indicates how “stretched” the curves are or, in other words, which order kinetics is involved in Eq. (12) for the structural parameter. Thus, in Paste 1 recovery curves jump up precipitously suggesting first-order kinetics, but, as we will see below, this is not always the case. After that we turn to the breakup curves, like those in Fig. 5, to determine the only remaining parameter—the breakup time Tb . Before we turn to a discussion of numerical solutions in vanecup geometry, it is worthwhile to compare simple estimates of relevant time scales. The time taken for the vane blades to rotate 2 90◦ is trot ∼/(2˝) = (1 − (R/R0 ) )/˙ app . Comparing this time to the recovery and breakup times from Eq. (9) we see that for all three pastes in the relevant ˙ app range between 0.01 and 30 s−1 we have trot  trec , tbreak . Hence, thixotropic effects are in fact important for these pastes since their structure fails to reach equilibrium between periodic hits by blades. Fig. 6 plots the same curves as in Figs. 4 and 5 and also adds the curves measured in vane-cup geometry. Here thick solid curves correspond to Couette axisymmetric geometry. Thin solid curves in Fig. 6 were calculated for transient model in 3D vane-cup geometry. Vane produces lower torques than Couette at all shear rates except the lowest, in which case it actually ends up slightly higher in the long-time limit after recovery. This is in agreement with 3D calculations. Dashed lines in Fig. 6 show results of calculations for the 3D equilibrium model as defined in Section 3, i.e., as if recovery and breakup were instantaneous. Steady Fluent solver was used to obtain these values. In Couette geometry transient and equilibrium solutions converge in the long-time limit and therefore the latter is not shown. For vane-cup, on the other hand, equilibrium (steady)

solution produces significantly lower torque than the vane-cup transient solution. The dashed lines in Fig. 6 correspond to the same final apparent shear rates as in the transient experiment and calculations: 3, 10 and 30 s−1 . We shall discuss the difference between the equilibrium and transient solutions in greater details below. Figs. 7 and 8 show the results of the same experiments as in Fig. 6, but for Pastes 2 and 3. Calculations also followed the same protocol: we first fitted Couette data and then used the same fitting parameters to calculate the torque in the vane-cup geometry. Considerable difference exists between the equilibrium flow curves for three pastes as well as their thixotropic parameters. Thus, Paste 3 is much more Newtonian (lower yield stress, higher high-shear viscosity) and also has lower n than the other two pastes. Paste 2 and especially Paste 3 have “stretched” recovery curves (hence, higher nr ) and stronger dependence of recovery time on shear rate (higher mr ). Figs. 9–14 show viscosity and structure inside the fluid volume calculated for the step-shear experiments in which apparent shear rate changes from 30 to 0.1 (recovery) and then from 0.1 to 10 s−1 (breakup). All these calculations were performed for Paste 1, but geometries, models and sometimes initial conditions varied. Fig. 9 deals with the effect of pre-shear on viscosity inside the gap in two different implementations of Couette geometry, 3D (Fluent) and 1D (iterative). Viscosity in the median cross-section is plotted against radial coordinate during different stages of the breakup experiment. These data correspond to the torque decay in Figs. 5 or 6. As the breakup progresses, the low-viscosity area at the cylinder wall penetrates deeper into the gap, but at the same time viscosity in the stagnant area at the cup wall increases. This demonstrates the effect of initial conditions in the thixotropic model, in this case incomplete recovery during 3000 s of shearing at 0.1 s−1 . Thus, if we extend the latter period to 30,000 s, viscosity in the stagnant area turns out much higher—recovery more complete. However, this has little effect on the torque, curves in Figs. 5 and 6

Fig. 9. Step-shear recovery curves calculated by CFD for Couette axisymmetric geometry (solid curves) and by iterative integration in 1D geometry, as described in Section 3.3 (dashed curves) with apparent shear rate changed as shown on the plot. ˙ min = 10−6 s−1 .

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

307

Fig. 10. Calculated contour plots for viscosity (in Pa s) marked by shades of grey around vane in 3D geometry. These plots correspond to the end of breakup stage (3000 s) at apparent shear rate 10 s−1 as defined by the experimental protocol. Radial line A and axial line B are defined in the text.

Fig. 11. Calculated contour plots for viscosity (a and b) and structural parameter (c) in horizontal cross-section of Fig. 10. Vane rotates counterclockwise.

308

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Fig. 12. Calculated viscosity plots along radial line A (a) and axial line B (b) as shown in Fig. 10 for the breakup step-shear test at different moments and apparent shear rates as shown. Dashed and solid curves correspond to initial condition  in = 2 and 10, respectively. Vertical line indicates side (on a) and top (on b) surface of the imaginary “cylinder” encompassing the vane. Origin is in the center of the vane, z-axis is directed up.

Fig. 13. Same as in Fig. 12 but for the breakup test which follows the recovery test of Fig. 12. Lines correspond to sequential moments of breakup as shown (in s).

Fig. 14. Calculated viscosity plots as in Figs. 12 and 13 for both, 3D vane-cup (dashed curves) and 2D-axisymmetric Couette (solid curves) geometries at the end of the recovery step (3000 s). Along with viscosity plots for the thixotropic model (thin curves), the curves for the equilibrium model (thick curves) are shown.

would have only slightly moved up at short times (less than by 10% at 1 s and less than by 5% at 10 s). In vane-cup geometry, as opposed to Couette, breakup and recovery processes never can be completely separated as each fluid element undergoes a periodic sequence of recovery-breakup events. Fig. 10 shows contour plots for viscosity around the vane calculated in 3D geometry for our equilibrium (a) and thixotropic (b) models at the end point of the “breakup” experiment in which shear rate changes from 0.1 to 10 s−1 . One can see that breakup turns out stronger around the vane on (a), while thixotropic recovery results in viscosity being higher between the blades on (b) which translates into higher torque in the thixotropic model (see Fig. 6 and below). One can also notice that breakup (the area of low viscosity) spreads farther from the vane in the thixotropic model (b), than in the equilibrium model (a), i.e., thixotropy makes breakup less concentrated in the area around the blades but also pushes it deeper into the sample. The same results are shown in Fig. 11(a) and (b) in the middle cross-section horizontal plane, while (c) shows the corresponding plot for the structural parameter in the thixotropic model on which one can better see the “wake” behind the vane blade. To have a better picture of recovery and breakup in vane-cup geometry, we defined two lines, A and B, shown in Fig. 10(b). Radial

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Fig. 15. Relative thixotropic torque increase, i.e., ratio of the torque calculated for thixotropic model to the torque calculated for time-independent equilibrium model with the same equilibrium viscosity, plotted against apparent shear rate for the three pastes.

line A lies in the median cross-section between the blades. Axial line B also lies in the middle between the blades being 7 mm distant from the axis, i.e., it lies partly on the surface of the imaginary cylinder which encompasses the vane and partly protrudes above it. Figs. 12–14 plot viscosity against radial and axial coordinate along lines A and B, respectively. Fig. 12 shows radial (a) and axial (b) plots at different moments of the recovery stage of the experiment on which apparent shear rate was changed from 30 to 0.1 s−1 . The lowest curves on both plots show viscosity profiles at the end of pre-shear, before apparent shear rate is changed to 0.1 and recovery starts. Recall that in our thixotropic model, once the shear rate changes, viscosity immediately changes as well, according to Eq. (5), even before structure starts to recover. Hence, Fig. 12 shows that at 0.2 s after the change of apparent shear rate to 0.1 viscosity had jumped up everywhere in the gap. In the next 3000 s during the actual recovery it increases further. Fig. 12 also shows the effect of the initial conditions. As was mentioned in the end of the previous section, our transient calculations were performed with the initial condition for the structural parameter  =  in (constant throughout the sample). Fig. 12 demonstrates that the actual value of  in is not really important since a memory of it is practically erased on the pre-shear stage. Thus Fig. 12 shows curves calculated with  in = 2 (dashed) and  in = 10 (solid). At the end of the 1000 s pre-shear the difference between the two solutions has already disappeared in the median plane (a) but it survives deep in the stagnant area above the vane (b) and eventually completely disappears by the end of the recovery. Torque curves in Figs. 4–6 coincide in both cases. Radial plot for viscosity on the breakup stage of the experiment is shown in Fig. 13. This plot may be compared to a similar plot for Couette geometry in Fig. 9 (the lower curves). The difference is quite obvious—in vane-cup geometry breakup spreads deeper into the sample while the stagnant area is less pronounced. Fig. 14 demonstrates the same conclusion more directly by plotting the results for the two geometries, Couette and vane-cup, and two models, equilibrium and thixotropic, on the same plot. Thin curves calculated at the end of the breakup (i.e., after 3000 s) are in fact the same as the ending curves in Fig. 9 (from the lower set) and Fig. 13 for Couette and vane-cup, respectively. Thick curves represent the equilibrium model. While viscosity profile for Couette

309

Fig. 16. Relative thixotropic torque increase (as in Fig. 15) for Paste 1 (Tr = 40 s, Tb = 15 s, and other parameters in Table 1, empty symbols) compared to calculations with only one of the parameters altered: either faster recovery (Tr = 4 s instead of 40 s, filled diamonds) or faster breakup (Tb = 5 s instead of 15 s, filled triangles).

geometry approaches its equilibrium limit (in fact, after additional 10 ks it practically coincides with it), the vane-cup viscosity profile does not develop any further, i.e., does not approach the profile for the equilibrium model. The fact that the equilibrium model fails to approximate the thixotropic model in vane-cup geometry has two important implications. The first implication is that the equilibrium model underestimates torque even for the constant-velocity problem which we have already seen in Figs. 6–8. The degree of an error associated with the equilibrium model depends on the properties of the fluid. Fig. 15 shows this more directly by plotting relative increase of torque due to thixotropy against apparent shear rate for all three of our pastes. Here the relative torque increase is defined as the torque calculated in the long-time limit (either after breakup or recovery—the result does not depend on the initial conditions) related to the torque calculated for the equilibrium model. One can see that for Paste 1 the torque increase due to thixotropy reaches 20% while for the other pastes it is maximized at about 10%. Exact value of the increase in each case depends on the competition between breakup and recovery. All other things being equal, faster recovery (smaller Tr ) and slower breakup (larger Tb ) both lead to stronger increase, as illustrated by Fig. 16. However, other parameters also have significant impact on the torque increase. For instance, Paste 1 has stronger torque increase in Fig. 15, as compared to Pastes 2 and 3, despite having slower recovery and faster breakup. This happens due to its steeper equilibrium curve, eq () ˙ (larger y and n) and stronger f multiplier (larger ˙ 0 ). These conclusions were verified by calculations which will not be discussed here in details. In general, there seems to be no simple rule to predict the torque increase without actually going through CFD calculations. Another implication of the equilibrium model failure is that the torque measured in the vane-cup geometry often turns out higher than could have been expected based on the measurement with Couette geometry and estimates based on this model. In Fig. 17 we plot the ratio of the torque calculated for vane-cup to the torque calculated for Couette for the same three pastes. Again, these results pertain to the constant-velocity problem. As is well known, equilibrium model predicts that vane-cup and Couette converge as yield stress approaches the yield point, but at high-shear rate the ratio decreases, i.e., vane-cup underestimates torque. Once we take thixotropy into account, the situation becomes more complex. First, for all pastes the ratio of computed vane-cup torque to Couette torque is higher, especially at high-shear rates. Second, for some

310

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

Fig. 17. Calculated ratio of the vane-cup torque to the Couette torque plotted against apparent shear rate for the three pastes. Empty symbols correspond to calculations for the equilibrium model; filled symbols—thixotropic model.

pastes (in this case, for Paste 1) it is actually higher than one, i.e., vane-cup gives higher torque than Couette. Overall, one can say that thixotropy in most cases reduces the difference between the two geometries. 5. Concluding remarks This paper dealt with two related problems: CFD modeling of mixing and Couette rheometry (with either coaxial cylinders or vane-cup setup). Mixing was not addressed directly in that we have not modeled any actual production mixers, but instead looked on a vane-cup rheometer as the simplest mixer. Both problems are well understood and can be easily solved for many non-Newtonian fluids, including yield-stress fluids, but thixotropic effects were never extensively analyzed, except for planar large-gap Couette with coaxial cylinders [7]. In this work we for the first time performed 3D computations of the fluid dynamics of vane-cup geometries for a thixotropic fluid. We want to reiterate that our goal was not to propose a new thixotropic model, but rather to demonstrate how any of such models (and there are quite a few of them) can be used for CFD calculations and why it is important, particularly for mixing. The model which we use is not entirely new—it may be thought of as a modification of structural models described in numerous reviews on thixotropy [1]. We had to make some modifications to make it suitable for our particular materials, toothpastes, but very similar results could have been obtained with many other models. The most important requirement in selecting the model was that it had to be suited for implementation in a commercial CFD code, such as Fluent. The most agreeable in this respect are the models relying on structural parameter(s) which may be treated as scalars obeying transport equations. Other important requirements for the model to be useful in industry were its flexibility and affordability of its implementation. For instance, toothpastes differ a lot in their equilibrium viscosity curves as well as transient behavior. Therefore, the model had to rely on multi-parameter equations in which the number of parameters was not particularly important. Furthermore, often the only available and reliable geometry for rheological characterization

of many pastes is large-gap Couette or vane-cup. Therefore, the model parameters had to be defined so that one could extract them from such measurements via fairly simple fitting procedure. In other words, our tasks were mostly computational, i.e., (i) how to extract quantitative data from practically affordable experiments (we used large-gap Couette to this end) and (ii) how to use these data to better model the flow of thixotropic fluids in mixers, of which a vane-cup rheometer itself may be thought of as the simplest example. The first task did not really require CFD (except, perhaps, to ascertain end effects) since we easily reduced it to integrals. The second was solved using a commercial CFD package, Fluent. The model we selected described quite well transient shear data for both Couette and vane-cup. This itself is not such a great achievement since other thixotropic models could have done that. Our main findings pertain to comparison between thixotropic and equilibrium models in general and between vane-cup and Couette. The first finding pertains to inadequacy of time-independent equilibrium models as compared to thixotropic models in describing mixers or rheometers. We found that thixotropy may significantly increase the torque on the shaft of a mixer (or vanecup rheometer for that matter). This happens because structure of a thixotropic fluid may not have enough time to reach equilibrium between moving blades. Consequently, the torque on the shaft of a mixer or rheometer will be underestimated if any of standard time-independent models are selected for the fluid. The degree of the underestimation depends on the properties of the fluid and is hard to predict without actually going through CFD calculations. Generally, the error reaches maximum at a certain angular velocity and was found as high as 20% for some pastes. The second finding pertains to the comparison between vanecup and Couette rheometers. For time-independent equilibrium models (i.e., neglecting thixotropy) calculating torque gives lower torque result for a vane-cup geometry than for an equivalent Couette geometry, with the relative difference increasing with shear rate. Since taking into account thixotropy increases torque for the vane-cup, it may either make this difference less pronounced or, for very thixotropic pastes, actually make torque for vane-cup higher than for Couette. Since practically all yield-stress materials are actually thixotropic, we expect these conclusions to be of a general nature. This may explain why vane-cup rheometer often performs better than could have been expected based on calculations typically performed for equilibrium models. There are some obvious limitations to our model, which may be addressed in the future. For instance, the model does not include viscoelasticity. The reason for this is primarily computational–the author is not aware of any software which would be capable of incorporating any viscoelastic model on a large enough scale to be applicable to industrial mixers, which is the ultimate purpose of this modeling. In Couette geometry one can easily generalize iterative procedures of Section 3.3 to include viscoelasticity, see Appendix B, but carrying this over to 3D seems an insurmountable task at the moment. Another limitation is that this model does not account for a totally gelled situation, i.e., it cannot be applied to the flow which starts from rest as reflected by the divergence of thixotropic relaxation times at low-shear rates. Nonetheless, the existing model is useful for modeling mixing and other industrial processes, which often are performed so that the gelled case is never approached. Appendix A. End effects Fig. 18 plots calculated correction factor for the constantvelocity problem in Couette geometry for equilibrium model for two pastes, Paste 1 and Paste 2. Thixotropy does not affect results in

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

311

which is similar to the horizontal cross-section in Fig. 2. Equilibrium model predicts the two estimates to converge in the low-shear limit to the same estimate (19), which is 1.19, and then drop to the noshaft estimate 1 + 2R/3L = 1.16. As we turn on thixotropy, the ratio increases roughly mirroring the overall torque increase (Fig. 15). Thus, the stronger is the impact of thixotropy on the torque, the stronger are end effects. Appendix B. Viscoelastic effects in 1D In simple shear instead of Eq. (11) we have for a Maxwell model: 

d + = f () ˙ , ˙ dt

(A.1)

where  is relaxation time. The solution (13) is generalized as = 2 + ( 1 − 2 )e−t/



− 2 Fig. 18. Calculated ratios of the torque in Couette 2D-axisymmetric model (T) to the torque in 1D model of infinite cylinders (T1D) and of T to the side contribution in the axisymmetric model (Tside), both plotted against apparent shear rate for Paste 1 and 2.

1−

eq (˙ 1 ) eq (˙ 2 )



tc g (e−t/tc − e−t/ )  − tc i

(A.2)

This equation predicts an undershoot minimum as the shear rate is reduced from high to low. For infinite concentric cylinders the iterative procedure similar to the one described in Section 3.3 can also be implemented. This is especially easy to do, if  is assumed constant. In this case instead of wall in Eqs. (16) and (18) we use the target stress, ¯ wall , to which wall approaches via relaxation, i.e.: ¯ wall (t) = 

d wall + wall . dt

(A.3)

Subsequently, Eq. (A.3) is solved together with (17) at each time step. These calculations, which will not be discussed here in details, allow somewhat better description of the initial segments of the curves for Couette geometry in Figs. 4–8 with  being an additional fitting parameter, but application of this model to vane-cup is problematic. References

Fig. 19. Calculated ratio of the torque in 3D to the torque in 2D vane geometries plotted against apparent shear rate for Pastes 1 and 2 in thixotropic (filled symbols) and equilibrium (empty symbols) models.

this case. The correction factor is properly defined as the ratio of the torque calculated in Couette 2D-axisymmetric model to the torque calculated in the 1D model. The Couette axisymmetric model also allows to separate torque contributions from different parts of the internal cylinder, for instance, from its side surface. Thus, we also plot of Fig. 18 the ratio of the complete torque to this side-surface contribution. In the low-shear limit both estimates approach (19), which for our particular system gives 1.19. As shear rate increases, the 2D-to-1D ratio first drops to 1.16, which is close to the noshaft estimate, 1 + 2R/3L, and then slightly increases at shear rate above 1 s−1 . Overall, we see that at least within 4% accuracy one can use a no-shaft estimate in the experimental shear rate range, while between 0.01 and 10 s−1 , the error does not exceed 2%. The sidesurface contribution actually drops at high-shear rates, suggesting that the observed increase of the torque over its 1D estimate at higher rates cannot be attributed to the side-surface flow, but rather has something to do with the flow near the edge of the cylinder. Similar plot for the vane-cup geometry is given in Fig. 19. In this case we are dealing with the ratio of the torque in the 3D model to the torque in a 2D vane model. The latter model is solved on a mesh

[1] H.A. Barnes, Thixotropy—a review, J. Non-Newtonian Fluid Mech. 70 (1997) 1–33. [2] K. Dullaert, J. Mewis, Thixotropy: build-up and breakdown curves during flow, J. Rheology 49 (2005) 1213–1230; K. Dullaert, J. Mewis, Stress jumps on weakly flocculated dispersions, J. Colloid Interface Sci. 287 (2005) 542–551. [3] H.A.B. Barnes, Q.D. Nguyen, Rotating vane rheometry—a review, J. nonNewtonian Fluid Mech. 98 (2001) 1–14. [4] Q.D. Nguyen, D.V. Boger, Yield stress measurement for concentrated suspension, J. Rheology 27 (4) (1983) 335–347. [5] H.A. Barnes, J.O. Carnali, The vane-in-cup as a novel rheometer geometry for shear-thinning and thixotropic materials, J. Rheology 36 (1990) 841–866. [6] S. Savarmand, M. Heniche, V. Bechard, F. Bertrand, P.J. Carreau, Analysis of the vane rheometer using 3D finite element simulation, J. Rheology 51 (2007) 161–177. [7] N. Roussel, R. Le Roy, P. Coussot, Thixotropy modeling at local and macroscopic scales, J. non-Newtonian Fluid Mech. 117 (2004) 85–95. [8] F. Moore, Rheology of ceramic slips and bodies, Trans. Br. Ceramic Soc. 58 (1959) 470–494. [9] D. Quemada, Rheological modeling of complex fluids. VI. Thixotropic and “thixoelastic” behavior. Start-up and stress relaxation, creep tests and hysteresis cycles, Eur. Phys. J. AP 5 (1999) 191–207. [10] C. Tiu, D.V. Boger, Complete rheological characterization of time-dependent food products, J. Texture Stud. 5 (1974) 329–338. [11] D. De Kee, R.K. Code, G. Turcotte, Flow properties of time-dependent foodstuffs, J. Rheology 27 (1983) 581–604. [12] E.A. Toorman, Modelling the thixotropic behavior of dense cohesive sediment suspensions, Rheol. Acta 36 (1997) 56–65. [13] M.R. Barkhadurov, C.I. Bronisz, C.W. Hirt, Three-dimensional thixotropic flow model. Published as report FSI-96-00-1 on www.flow3d.com (1996). [14] T.C. Papanastasiou, Flow of materials with yield, J. Rheology 31 (1987) 385–404. [15] A.A. Potanin, R. De Rooij, D. Van den Ende, J. Mellema, Microrheological modeling of weakly aggregated dispersions, J. Chem. Phys. 102 (1995) 5845–5853.

312

A. Potanin / J. Non-Newtonian Fluid Mech. 165 (2010) 299–312

[16] A. Mujumdar, A.N. Beris, A.B. Metzner, Transient phenomena in thixotropic systems, J. non-Newtonian Fluid Mech. 102 (2002) 157–178. [17] K. Dullaert, J. Mewis, A structural kinetic model for thixotropy, J. nonNewtonian Fluid Mech. 139 (2006) 21–30. [18] C. Macosco, Rheology. Principles, Measurements and Applications, Wiley, 1994.

[19] C. Baravian, D. Quemada, A. Parker, Modeling thixotropy using a novel structural kinetics approach: basis and application to solution of iota carrageenan, J. Texture Stud. 27 (1996) 371–390. [20] O. Wein, M. Vecer, J. Havlica, End effects in rotational viscometry. I. Noslip shear–thinning samples in the Z40 DIN sensor, Rheol. Acta 46 (2007) 765–772.