Thixotropic yield-stress fluids (TYSFs) are very common and include emulsions [1], pastes [2], personal care products [3], colloidal suspensions [4,5], and foodstuffs [6–8]. TYSFs usually have weakly aggregated internal structures that break under flow and can be restored by the Brownian motion or interactions involving relative motion of particles [9]. Structure breakup and reformation, also referred to as shear rejuvenation and aging, lead to the time-dependent and reversible rheological behaviors of TYSFs. The rheological responses of TYSFs often involve both viscoelastic and thixotropic components. Just as accurate descriptions of viscoelasticity usually require multiple relaxation times—a reflection of internal structural evolution at multiple length scales—quantitative predictions of thixotropic response require multiple modes to capture structural kinetics. The existence of multiple thixotropic timescales can be revealed through non-monotonic time-dependence of rheological responses in multiple-step shear-rate-jump [10] or stress-jump [4] tests. The multiple viscoelastic relaxation times of fluids can be identified by fitting an empirical model to experimental data from small amplitude oscillatory shear tests at different frequencies. However, the small deformation used to measure the relaxation times does not break the internal structures and so does not characterize the thixotropic character of the fluid. To measure multiple thixotropic relaxation times, one must turn to transient shear tests such as shear-rate-step tests. Of course, a fluid can have both viscoelastic and thixotropic characters; our focus here is on the latter.

We previously [10,11] proposed a technique, which we call the “multi-lambda approach” (MLA), to generalize thixotropic models that would otherwise involve only a single structure parameter λ. This approach introduces multiple structure parameters which we denote as λi. The λis relaxation rates combine to produce a stretched-exponential relaxation. The MLA requires only a single additional model parameter in addition to those involved in the corresponding single lambda model; the multilambda approach improves predictions of transient thixotropic responses [10,11].

In this letter, we apply the MLA to the model of Coussot et al. [4] to predict the nonmonotonic time evolution of shear rate in multiple-step stress-jump tests. Coussot et al. [4] suggested introducing multiple thixotropic timescales to predict this phenomenon; however, no such introduction has been available or reported to date. Although other data sets [12,13] might be equally well suited to application of the MLA, we here select the Coussot et al. [4] data set because it highlights the interesting nonmonotonic viscosity bifurcation. Such behavior is a hallmark of the phenomenology that the MLA is well-suited to model quantitatively.

Coussot et al. [4] studied the rheological behaviors of a Bentonite suspension (solid volume fraction of 4.5%) in a set of stress-jump experiments. The test protocol includes (1) preshearing at a relatively large stress producing a large steady shear rate (400s1), (2) resting for a duration of tw, and (3) shearing again at different lower stress levels than the preshear stress. Figure 1(a) shows the results of the experiments with tw=0. After the sudden decrease of the shear stress, the viscosity gradually builds up and therefore the shear rate decreases. The viscosity bifurcates at about 6 Pa. That is, for higher stresses than 6 Pa, the material flows continuously with the shear rate slowly decreasing. For stresses lower than this, the shear rate eventually falls rapidly below 104s1; i.e., flow effectively ceases. The rheological responses are drastically different if there is a short delay time tw between the preshear and the final step, as shown in Fig. 1(b) for tw=20s. In this case, the viscosity bifurcates at about 12 Pa so that for stresses lower than this value, flow eventually stops. Interestingly, for stresses within the range (6–12 Pa), the shear rate first increases and then rapidly collapses, indicating that at a constant stress structural breakage and growth can coexist—the former dominates the short-term relaxation and the latter dominates the long-term evolution. Note that here we describe the nonmonotonic evolution of the shear rate as a thixotropic response. It is very unlikely to be predominantly elastic because both the time scales (around 10 s) and the strain scales (around 100 strain units) are much larger than those of a typical elastic response (t < 0.1 s and γ < 0.01).

Coussot et al. [4] proposed a simple constitutive model for viscosity-bifurcating TYSFs that qualitatively captures most of their experimental findings except for the nonmonotonic time evolution of shear rates shown in Fig. 1(b). Their model involves the stress expression,

σ=η0f(λ)γ˙,
(1)

where σ denotes the shear stress, η0 the asymptotic value of the viscosity, γ˙ the shear rate, λ the structure parameter, and f(λ) an unknown function of λ. f(λ) satisfies two constraints: (1) it is an increasing function of λ, and (2) f(0)=1, where λ is zero for a completely broken structure. Coussot et al. [4] found that f(λ)=exp(λ) qualitatively agrees with their experiments. Their expression for the evolution of λ is

dλdt=1T0αλγ˙,
(2)

where T0 is the characteristic timescale of aging and α is a model parameter. In this model, the value of λ is unbounded. When γ˙=0, λ continuously builds up to infinity. Equations (1) and (2) can be written in dimensionless forms by defining the dimensionless stress Σ=σαT0/η0, the dimensionless shear rate Γ=αT0γ˙, and the dimensionless time Θ=t/T0,

Σ=exp(λ)Γ,
(3a)
dλdΘ=1λΓ.
(3b)

At steady state, the flow curve follows

Σ=Γexp(Γ1).
(4)

The flow curve has a negative slope for Γ<1. In a constant stress test at Σ=Σ0, the time evolution of λ follows

dλdΘ=1λexp(λ)Σ0,
(5)

where λ always evolves monotonically—either toward the steady state value at Σ0 or toward infinity. Regardless of the initial value of λ, a nonmonotonic evolution of λ requires dλ/dΘ to have nonunique values at certain λ values, which is impossible for Eq. (5). At Σ=Σ0, Γ(Θ)=Σ0/exp[λ(Θ)]. Γ(Θ), therefore, also evolves monotonically, and the model cannot predict the nonmonotonic behavior shown in Fig. 1(b).

To predict the nonmonotonic phenomenon, we generalize this simple model using the MLA [10,11]. We decompose λ into a group of substructure parameters λi according to

λ=i=1NCiλi,
(6)

where Ci is the weight coefficient of λi and N is the number of thixotropic modes. The evolution equation for each λi is

dλidΘ=Di(1λiΓ).
(7)

Here, Di is a dimensionless rate coefficient of λi. Note that N, Ci, and Di are not model parameters, but they are obtained from a stretched exponential distribution according to

exp(sβ)i=1NCi(β)exp(Dis),
(8)

where s is a dummy time variable and β is the stretching exponent, which is the only new model parameter we introduce. Equation (8) states that a stretched exponential function can be approximated by a linear combination of simple exponentials. This decomposition is, by itself, a numerical problem and is independent of other equations in the model. The choice of a stretched-exponential distribution is based on a previous experimental finding that the stress evolution after a single-step-rate jump follows a stretched-exponential function [10]. Note that in this new model, in single-step-rate tests, λ(t) follows a stretched exponential function, but σ(t) does not because σ(t)exp[λ(t)]γ˙. Nevertheless, when the change of λ(t) over time is small, e.g., λ(0)/λ()(0.1,10), σ(t) can be closely approximated by a stretched exponential.

Equation (3a) and Eqs. (6)–(8) specify the dimensionless form of this new model. When β=1, N=Ci=Di=1, and the new model reduces to the original form of Eq. (3). When β<1, one can use arbitrarily many modes for the approximation described by Eq. (8); however, the smaller β is, the more modes are needed for an accurate representation of the stretched exponential. For β=0.5, 10 modes can give a very close approximation [11]. The method for calculating Ci and Di for a given value of β can be found in our previous works [10,11]. We reiterate that even with the selection of 10 modes, the MLA still only adds a single parameter, namely, β, to the set that must be specified.

We show the predictions of the new model in Figs. 1(c) and 1(d). Here, we focus on a parsimonious formulation and hence make only qualitative comparisons between this new model and the experiments. We arbitrarily set β=0.5 and use 10 modes for the thixotropic kinetics. Figure 1(c) shows the predictions of the dimensionless shear rate Γ(Θ) after Σ instantly jumps from a high value to several lower levels. The impact of multiple structure parameters is not obvious in this case because Θw=0 and, therefore, all λi values evolve upward. The influence of multiple structure parameters can be clearly seen if a short rest period is introduced between the preshear and the subsequent lower stress. During rest, the λis with large Di values rapidly build up, and after restart of the stress they break down, causing the initial rise of Γ. Those λis with small Di values, however, do not build up much during resting. They have longer memories of the flow history. At stress levels lower than the preshear stress, they continue to slowly build up and cause the long-term decrease of Γ. The predictions of this simple model agree nicely with experiments, as shown in Fig. 1.

Summary. In this letter, we introduced multiple modes of thixotropic kinetics to the model of Coussot et al. [4] for a viscosity-bifurcating yield-stress fluid to predict the nonmonotonic time evolution of shear rates in multiple stress-jump tests that their original model cannot capture. Our modification requires only one additional model parameter and allows the original model to now qualitatively capture the phenomenon. This letter, together with our previous studies [10,11], strongly demonstrates the utility of including multiple structural kinetics in thixotropic models. This inclusion is especially important for complex shear histories that extend the applicability of TYSF models beyond tests such as step jumps in shear rate or stress. We note that this extension is phenomenological and provides no microstructural explanation for the multiple thixotropic time scales whose origins would require a microscopic theory to address.

The authors thank Procter & Gamble Company for funding the work of Y.W., and R.G.L. acknowledges funding by the National Science Foundation (NSF), under Grant No. CDS&E-1602183. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors and do not necessarily reflect the views of the NSF.

1.
Paredes
,
J.
,
N.
Shahidzadeh-Bonn
, and
D.
Bonn
, “
Shear banding in thixotropic and normal emulsions
,”
J. Phys. Condens. Matter
23
(
8
),
284116
(
2011
).
2.
Durairaj
,
R.
,
N. N.
Ekere
, and
B.
Salam
, “
Thixotropy flow behaviour of solder and conductive adhesive pastes
,”
J. Mater. Sci. Mater. Electron.
15
,
677
683
(
2004
).
3.
Kropf
,
C.
, and
U.
Brueninghaus
, “
Thixotropic oral and dental care agents
,”
Google Patent US
20030175217 A1 (
2001
).
4.
Coussot
,
P.
,
Q. D.
Nguyen
,
H. T.
Huynh
, and
D.
Bonn
, “
Viscosity bifurcation in thixotropic, yielding fluids
,”
J. Rheol.
46
,
573
589
(
2002
).
5.
Fall
,
A.
,
J.
Paredes
, and
D.
Bonn
, “
Yielding and shear banding in soft glassy materials
,”
Phys. Rev. Lett.
105
,
225502
(
2010
).
6.
Augusto
,
P. E. D.
,
V.
Falguera
,
M.
Cristianini
, and
A.
Ibarz
, “
Rheological behavior of tomato juice: Steady-state shear and time-dependent modeling
,”
Food Bioprocess Technol.
5
,
1715
1723
(
2012
).
7.
Razavi
,
S. M. A.
, and
H.
Karazhiyan
, “
Flow properties and thixotropy of selected hydrocolloids: Experimental and modeling studies
,”
Food Hydrocolloids
23
,
908
912
(
2009
).
8.
Yousefi
,
A. R.
, and
S. M. A.
Razavi
, “
Steady shear flow behavior and thixotropy of wheat starch gel: Impact of chemical modification, concentration and saliva addition
,”
J. Food Process Eng.
39
,
31
43
(
2016
).
9.
Vermant
,
J.
, and
M. J.
Solomon
, “
Flow-induced structure in colloidal suspensions
,”
J. Phys. Condens. Matter
17
,
R187
R216
(
2005
).
10.
Wei
,
Y. F.
,
M. J.
Solomon
, and
R. G.
Larson
, “
A multimode structural kinetics constitutive equation for the transient rheology of thixotropic elasto-viscoplastic fluids
,”
J. Rheol.
62
,
321
342
(
2018
).
11.
Wei
,
Y.
,
M. J.
Solomon
, and
R. G.
Larson
, “
Quantitative nonlinear thixotropic model with stretched exponential response in transient shear flows
,”
J. Rheol.
60
,
1301
1315
(
2016
).
12.
Armstrong
,
M. J.
,
A. N.
Beris
,
S. A.
Rogers
, and
N. J.
Wagner
, “
Dynamic shear rheology of a thixotropic suspension: Comparison of an improved structure-based model with large amplitude oscillatory shear experiments
,”
J. Rheol.
60
,
433
450
(
2016
).
13.
Dullaert
,
K.
, and
J.
Mewis
, “
Thixotropy: Build-up and breakdown curves during flow
,”
J. Rheol.
49
,
1213
1230
(
2005
).