Letter to the Editor: Modeling the nonmonotonic time-dependence of viscosity bifurcation in thixotropic yield-stress fluids

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 $λi′s$ 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 $(≈400s−1)$⁠, (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 $10−4s−1$⁠; 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,

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

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$⁠,

At steady state, the flow curve follows

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

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

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

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

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 $λi′s$ with large $Di$ values rapidly build up, and after restart of the stress they break down, causing the initial rise of $Γ$⁠. Those $λi′s$ 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.