To predict the complex transient rheology of thixotropic elasto-viscoplastic (TEVP) fluids, we generalize a previous scalar thixotropic-only “multilambda” (ML) model [Wei *et al*., J. Rheol. **60**(6), 1301–1315 (2016)] and combine it with the isotropic kinematic hardening (IKH) model [C. J. Dimitriou and G. H. Mckinley, Soft Matter **10**(35), 6619–6644 (2014)]. This new constitutive equation, which we call ML-IKH model, has the following features: (1) Multiple thixotropic structure parameters that collectively exhibit a stretch-exponential thixotropic relaxation in step tests; (2) nonlinear thixotropic kinetic equations in both the shear rate or stress and the structure parameters; (3) incorporation of the Armstrong-Frederick kinematic hardening rule [C. O. Frederick and P. Armstrong, Mater. High Temp. **24**(1), 1–26 (2007)] for the evolution of yield stress; and (4) viscoelasticity. We evaluate this 12-parameter model, discuss its four key features, and compare its predictions with those of the ML, IKH, and modified Delaware thixotropic models [Armstrong *et al*., J. Rheol. **60**(3), 433–450 (2016)] for two sets of experimental data [Wei *et al*., J. Rheol. **60**(6), 1301–1315 (2016); Armstrong *et al*., J. Rheol. **60**(3), 433–450 (2016)] of a TEVP fumed silica suspension. The shear-rate histories include steady state, step shear rate, step stress, intermittent shear, flow reversal, and large amplitude oscillatory shear (LAOS). We show that in step tests the thixotropic and viscoelastic evolutions are dominant, while intermittent shear tests the multiple thixotropic timescales, and flow reversal tests the viscoelastic and plastic evolutions. The rheological responses in LAOS tests are more complex and involve all aspects of TEVP rheology. The four features quantitatively capture different aspects of TEVP rheology. We also provide a tensorial formulation of the ML-IKH model that is frame-invariant, obeys the second law of thermodynamics, and can reproduce the predictions of the scalar version.

## I. INTRODUCTION

Many complex fluids exhibit flow-sensitive rheological properties due to flow-induced changes of internal structures [1–3]. Complex fluids with weakly aggregated internal structures often display thixotropy, defined as a flow-induced reversible time-dependent shear-thinning (or sometimes thickening) viscous response [4–6]. The shear thinning is produced by deformation and breakage of internal structures under flow. Because of Brownian motion or interactions that induce relative motion of particles [4], the microstructure can rebuild as flow decreases or stops, and, in some cases, form a space-filling structure network. In these systems, the deformation of the structure may store and release elastic energy before breaking, which gives rise to viscoelasticity as well as thixotropy. The system also exhibits a yield stress if the space-filling structural networks can withstand a modest stress without breaking or continuously deforming. Most thixotropic fluids also simultaneously show both viscoelasticity and yield stress [7,8], making it difficult to develop accurate constitutive equations for them. Such a complex rheological response is often referred to as “thixotropic elasto-viscoplastic” (TEVP) [8].

Thixotropic fluids are very common and include blood [9], paints [10], adhesives [11], concrete [12], personal care products [13,14], waxy crude oil [15], food materials [16–18], and weakly flocculated colloidal suspensions [19–23]. The rheology of these fluids plays a key role in the design of processing flows, the development of products with a tailored thixotropic response, and the general study of microstructure-rheology relationships in soft matter.

A thixotropic constitutive model usually contains both an equation for the stress tensor and an auxiliary equation for the time evolution of one or more (usually scalar) structure parameters. For example, the models in [15,24–28] introduce a structure parameter *λ* that varies between zero and one, representing, respectively, a completely broken-down state and a fully structured one. (Some authors [29,30] use fluidity as the structure parameter, which has an inverse relationship to the size or strength of the structure.)

Most models use empirical relationships to link the structure parameter with bulk rheological properties. For example, the assumption of a linear dependence of viscosity and yield stress on *λ* is common. Others assume a specific microstructural physical picture and define the structure parameter and rheological properties based on physical modeling. For example, the relationship between viscosity and the effective floc volume fraction has been adopted in [31,32].

As discussed in the recent review [8], TEVP materials show complex rheology due to the coexistence of thixotropy, viscoelasticity, and plasticity. They show distinct rheological behaviors in test conditions with different amplitudes of strain, strain rate, and stress; and different timescales of thixotropic, viscoelastic, and plastic evolution. Previous thixotropic models, nevertheless, have focused on only limited test conditions. Although they achieve success in predicting certain tests, they usually lack the applicability in capturing the complex TEVP rheology. For example, our previous thixotropic-only multilambda (ML) model [33] successfully captures the nonlinear and stretched-exponential thixotropic evolution, but it lacks the component of viscoelasticity and therefore gives poor predictions in transient tests when the shear rate is small. The isotropic kinematic hardening (IKH) model shows good agreement with the rheology of waxy crude oils in several hysteresis tests, shear start-up tests, and large amplitude oscillatory shear (LAOS) tests. However, a comprehensive examination of the IKH model in other rheological experiments is lacking. The modified Delaware thixotropic (MDT) model was examined in a range of step-shear-rate tests, flow reversal tests, and LAOS tests. It successfully captures several important features of TEVP rheology including elastic stress contribution in steady state, thixotropic buildup and breakdown, and viscoelasticity in LAOS tests. But it lacks the feature of multiple thixotropic time scales, which is important to accurately predict the thixotropic response; it also shows poor predictions in flow reversal tests and some LAOS tests.

Still lacking, therefore, is a comprehensive model that can give quantitative predictions over wide-ranging test conditions. Developing such a robust model would be important in applications that may involve all the above-mentioned or more complex flow conditions. Another limitation of the available thixotropic models is that they are written mainly in scalar forms and are therefore directly applicable only to simple shear flows. The relatively few [15,34–36] that have been written in tensorial form have for the most part not been tested against experimental results. Thus, a well-tested, quantitative tensorial constitutive model for thixotropic fluids is currently lacking.

This work takes some steps to fill this gap. In Sec. II, we review and discuss four critical attributes that TEVP models should possess. Then in Sec. III, we introduce the ML-IKH model for TEVP fluids. We present the model in both a scalar and a tensorial form, the latter of which is in Eulerian form, frame-invariant, and compliant with the second law of thermodynamics. We outline the procedure for model parameterization in Sec. IV. And in Sec. V, we test the proposed model, in comparison with the ML, IKH, and MDT models, over a wide range of simple shearing flows, including steady state shear, unidirectional step tests, intermittent step tests, flow reversal tests, and LAOS tests. (We leave tests under nonshearing kinematics to the future.) We discuss the model components that are critical to successful prediction of those rheological tests. In Sec. VI, we discuss advantages of the tensorial model. We finally summarize the results and discussion in Sec. VII.

## II. THEORETICAL OVERVIEW

### A. Predicting ideal thixotropy

In structural kinetics models of thixotropy, a structure parameter, usually denoted as *λ*, is introduced to indicate the instantaneous degree of “structure.” The constitutive relationship between stress and the rate of deformation is parameterized by *λ* whose rate of change is assumed to be a function of the instantaneous flow condition and *λ* itself. While almost all thixotropic fluids exhibit viscoelasticity and yield stress, their rheological behavior is close to that of purely viscous thixotropic fluids when the shear rate is sufficiently large.

As discussed in our previous work [33], a simple and general model of ideal thixotropy is

Here, *σ* denotes the shear stress, $\gamma \u0307$ the shear rate, *α* the order of dependence of stress on shear rate, *K*(*λ*) the coefficient which is an increasing function of *λ*. *k*_{1}, *k*_{2}, and *k*_{3} are, respectively, the rate coefficients of the break-down term, the shear-induced build-up term, and the Brownian build-up term; *a* and *b* are the orders of the dependence of the first two of these terms on the flow parameter *ϕ*; and *n* and *m* are the orders of the dependence of *λ* in the break-down term and build-up terms. The “flow parameter,” *ϕ*, can be shear rate, shear stress, or a combination of both. A summary of representative kinetic equations can be found in the review [4].

Most thixotropic models assume linear relationships between *K*(*λ*) and *λ*, *σ* and $\gamma \u0307$, d*λ*/d*t* and *λ*. Such linear thixotropic models predict that *λ*, and therefore *σ*, respond exponentially to step changes in *ϕ* (assumed to be $\gamma \u0307$ in most cases). However, experimental results [33,37,38] show that thixotropic relaxation often follows a stretched exponential function. In addition, monolambda thixotropic models in the form of Eqs. (1) and (2) predict that *λ*, if initially out of equilibrium, always relaxes monotonically to its steady-state value when the shear rate or stress is constant, regardless of any prior shear history. Consequently, the stress, viscosity, and/or shear rate evolve monotonically. However, nonmonotonic thixotropic relaxation has been observed at constant shear rate, if this constant shear rate follows a previous multistep flow history. Such nonmonotonicity indicates the coexistence of multiple time scales of thixotropic relaxation that single-lambda thixotropic models fail to predict [33].

We previously developed the ML approach [33] to incorporate a spectrum of structure parameters into single-lambda thixotropic models at the cost of only one additional model parameter. This approach starts by writing *λ* as a linear combination of a spectrum of “substructure-parameters,” *λ _{i}*

where *C _{i}* is the coefficient

*λ*and

_{i}*N*is the number of thixotropic relaxation modes. Equation (2) is replaced by the kinetic equation of each

*λ*

_{i}According to Eq. (4), all *λ _{i}* have the same steady-state value as

*λ*in Eq. (2) at the same shear rate or stress; however, their time derivatives differ and are set by the dimensionless prefactor

*D*. Equations (3) and (4) imply that the thixotropic evolution can be divided into multiple independent modes with different relaxation times. Equations (3) and (4) reduce to Eq. (2) if only a single mode (

_{i}*N*= 1) is included and both ${Ci}$ and ${Di}$ have only a single element of unity.

To account for the stretched exponential thixotropic evolution in single step tests, we assume that the thixotropic relaxation times follow a stretched exponential distribution—the ${Ci}$ and ${Di}$ are chosen to satisfy the following relationship:

where *β* is the stretching exponent, 0 < *β < *1. Equation (5) shows the decomposition of a stretched exponential function into a linear combination of simple exponentials. This decomposition is, by itself, a numerical problem and is independent of other equations in this model. It requires the input of *β* and outputs the values of *N*, ${Ci}$, and ${Di}$. Thus, *β* is the only additional model parameter accounting for the stretched exponential thixotropic evolution. For a given *β*, *N*, ${Ci}$, and ${Di}$ are determined through the following procedure. First, the user specifies *N* based on the requirement of the approximation. *N* sets the resolution of the decomposition. To achieve a certain level of approximation, the smaller *β* is, the larger *N* needs to be. *β = *1 indicates a simple exponential (*N* = *C* = *D* = 1). Second, ${Di}$ is generated containing *N* modes with a range of relaxation times. Last, the optimal ${Ci}$ is calculated to best approximate the stretched exponential. Details regarding this decomposition (the mathematical background, the algorithm for calculating $Ci$ and $Di$, the influence of the choice of *N* and $Di$) can be found in our previous work [33]. An improved subroutine to perform this decomposition is provided in the supplementary material [39].

### B. Incorporation of viscoelasticity into thixotropic models

Many thixotropic fluids are simultaneously viscoelastic. Viscoelasticity manifests itself in step-shear-rate tests where the stress may first exhibits an instantaneous jump, followed by a relatively fast (viscoelastic) relaxation that subsequently transitions to a slower (thixotropic) one in the opposite direction. For example, upon a stepwise increase of shear rate, stress may jump up instantly, increases toward a maximum and then gradually decays to a plateau higher than its initial steady-state value. The contribution of elastic stress to the total stress can be directly measured in flow cessation tests where flow suddenly stops after a period of preshear [40]. The residual stress that does not instantly vanish when shear ceases is the elastic stress.

Thixotropic models need to incorporate viscoelasticity to capture such behaviors. Figure 1 is a summary of representative approaches. Figure 1(a) is derived from the Bingham model, with the viscosity written as a function of *λ*, and the yield stress replaced with an elastic stress (*σ*_{e}) that is a function of both *λ* and $\gamma \xafe$. $\gamma \xafe$ is a strainlike structure parameter. An additional kinetic equation is needed to describe the evolution of $\gamma \xafe$. Examples of constitutive equations of the type of Fig. 1(a) are the models of Mujumdar *et al*. [26], and Dullaert and Mewis [27].

Figure 1(b) is derived from the Maxwell model by assuming a *λ*-dependent viscosity, $\eta (\lambda )$, and/or modulus, $G(\lambda )$, each of which contributes to the strain, $\gamma $, which is additively divided into the viscous and elastic strains, $\gamma v$ and $\gamma e$, respectively. Examples of models of this form were proposed by Labanda and Llorens [41], de Souza Mendes [42], and Dimitriou and McKinley [15]. Figure 1(c) is derived from the Oldroyd-B (or Jeffreys) model, in that it contains a parallel viscous Newtonian stress that is structure independent, along with the structure-dependent viscoelastic Maxwell component. Examples are the models of de Souza Mendes and Thompson [43], Blackwell and Ewoldt [44], and the model in this work. Figure 1(d) is a hybrid of the Bingham- and Maxwell-type models which assumes an additive decomposition of both stress and strain. The elastic stress is a function of the elastic strain; and the viscous stress is a function of the viscous shear rate, with both also functions of *λ*. Examples are the models of Armstrong *et al*. [24] and Mwasame *et al*. [31].

In this work, we propose a model in the form of Fig. 1(c), which naturally captures the following rheological features resulting from viscoelasticity: (1) An instantaneous stress jump when the shear rate undergoes a stepwise change, which arises from the medium viscosity *η*_{m}, (2) a continuous but fast initial elastic relaxation resulting from the spring-dashpot element containing *G* and *η*, and (3) a slower thixotropic relaxation in the direction opposite to that from the elastic relaxation in step-shear-rate tests, which results from the *λ*-dependence of *G* and *η*. There is also a gradual, but fast, decaying residual stress after cessation of shear, resulting from the spring-dashpot element.

### C. Kinematic hardening

“Kinematic hardening” during yielding is commonly invoked to account for the Bauschinger effect [49] that many materials, especially metals, exhibit [45–49]. The Bauschinger effect refers to an anisotropic yield stress resulting from plastic deformation. Kinematic hardening models describe this phenomenon through a changing yield surface, the center of which moves as plastic deformation occurs, such that the material becomes harder to deform (exhibits a larger yield stress) in the direction of past deformation, but easier to deform in the reverse direction.

While extensively employed in the plasticity literature, kinematic hardening has only been recently introduced into thixotropic modeling, by Dimitriou and McKinley [15]. Their IKH model incorporates a modified Armstrong-Frederick kinematic hardening rule [49] into a Maxwell-type thixotropic model. This kinematic hardening component was found to be critical to the ability of the IKH model to predict the dynamic rheology of a waxy crude oil in LAOS tests [15].

The Armstrong-Frederick kinematic hardening rule is widely used in the limit of small strains to model the Bauschinger effect. Incorporating the Armstrong-Frederick kinematic hardening rule into the Bingham model gives Eqs. (6)–(9). First, a back stress ($\sigma back$) is introduced. Here, *A* is a structure parameter often called the “back strain”; and *k*_{h} is a constant coefficient

Equation (7) describes the evolution of *A*

where *q* is a material constant. Equation (7) implies that *A* is determined solely by the strain histories. $\sigma back$ is equivalent to an elastic stress and *A* an elastic strain when *q* = 0. Next, Eq. (8) defines the effective stress ($\sigma eff$) and Eq. (9) gives the relationship between $\gamma \u0307$ and $\sigma eff$.

where *k*_{y} is the yield stress in the initial undeformed state where *A* = 0 and $\eta $ is the plastic viscosity.

The Armstrong-Frederick kinematic hardening rule introduces a varying yield surface. Initially in the undeformed state where *A* = 0, the yield surface in stress space is $[\u2212ky,ky]$. As strain occurs and *A* evolves toward the direction of strain, the yield surface moves to $[\u2212ky+khA,ky+khA]$. The yield stress in the direction of deformation therefore increases and that in the opposite direction decreases. This model reduces to the Bingham model if *k*_{h} = 0. Equation (9) can be further generalized, giving Eq. (10) [50].

The above equation reduces to the original Armstrong-Frederick equation when *m* = 1. As reported by Dimitriou and McKinley [15], this modification is important in predicting the dynamic rheology of waxy crude oils (*m* = 0.25). We find *m* = 1 is in good agreement with the rheology of the fumed silica suspension.

In this work, we introduce the Armstrong-Frederick kinematic hardening rule into the proposed model and test it against the rheological response of a model thixotropic fluid [23], a weakly aggregated fumed silica suspension. We show that the kinematic hardening rule is critical to predicting the rheology of this fluid in flow reversal tests, LAOS tests, and other complex flow histories involving a changing deformation direction.

## III. PROPOSED MODEL

Section II discussed four important aspects of modeling the rheology of TEVP fluids, namely: (1) Multiple thixotropic structure parameters, (2) nonlinear thixotropic evolution equations, (3) viscoelastic model components, and (4) a yield stress with kinematic hardening. In this section, we introduce the ML-IKH model which includes these four aspects. (In comparison, the ML model includes only aspects 1 and 2; the IKH model contains only aspects 3 and 4; and the MDT model only has aspect 3.) This model generalizes our previous ML “ideal thixotropic model” [33] by amalgamating it with the IKH model [15], and making some additional modifications, as discussed shortly. We propose the model in both scalar and tensorial forms.

### A. Scalar form

The total stress is composed of contributions from the medium (*σ*_{m}) and from internal structures (*σ*_{s})

In this work, we consider only thixotropic suspensions with a Newtonian continuous phase. Assuming that the system has a constant viscosity of *η*_{m} when all internal structures are broken, we write *σ*_{m} as follows:

As a common assumption in elasto-viscoplastic models [15,24,51], the total shear strain is additively decomposed into the elastic strain (*γ*_{e}) and plastic strain (*γ*_{p})

Here, *γ*_{e} represents the elastic deformation of the internal structures and *γ*_{p} represents the plastic distortion caused by, for example, the relative movement of internal structures with respect to the medium. $\gamma \u0307e$ and $\gamma \u0307p$ are, respectively, the elastic and plastic strain rate. We assume linear elasticity and a constant modulus, taking *σ*_{s} is proportional to *γ*_{e}

Next, we introduce the equations determining $\gamma \u0307p$. $\gamma \u0307p$ can be written as the product of its sign (*n*_{p}) and its absolute value ($|\gamma \u0307p|$)

Here, *n*_{p} indicates the direction of the plastic deformation while $|\gamma \u0307p|$ is the magnitude of its rate. Both *n*_{p} and $|\gamma \u0307p|$ are determined by *σ*_{eff}. Slightly different from Eq. (8), here *σ*_{eff} is defined as

As is common in plasticity theory [52], we assume that the plastic deformation is codirectional with the effective stress, i.e.,

$|\gamma \u0307p|$ is determined as follows:

where *η*_{thi} is the plastic viscosity.

Equations (11)–(18) can be written in a compact form as follows:

where $\theta =\eta thi/G$ is the viscoelastic relaxation time in the fully structured state. Notice that when $ky=0$, $kh=0$, and $\lambda $ is removed from this model (by setting $\lambda =1$), this model reduces to the Jeffrey model; when $\theta =0$, $kh=0$, and $\lambda $ is absent, this model reduces to the Bingham model.

We next introduce the evolution equations for $A$ and $\lambda $ to complete this model. Equation (7) is modified to relate the evolution of back strain to the plastic strain rate

Following the ML approach, Eqs. (3)–(5) are incorporated to describe the evolution of $\lambda $.

It has been shown [24,27,33,53] that the thixotropic build-up rate is linear in $\lambda $ (*m* = 1) if a linear relationship between $\lambda $ and viscosity or yield stress is adopted. Previous studies [27,33] in addition show that the flow-induced aggregation rate has a square-root dependence (*b* = 0.5) on the shear rate or stress. This relationship is based on the theoretical study of van de Ven and Mason [54] concerning the effect of shear flow on the aggregation kinetics of colloidal dispersions at low Peclet numbers. The build-up terms with *m* = 1 and *b* = 0.5 were first adopted by Dullaert and Mewis [27] and used by Wei *et al*. [33].

The parameters *a* and *n* in the break-down term of Eq. (4) are two adjusting model parameters. The nonlinearity of *λ _{i}* in the breakdown term captures the positive correlation between the thixotropic relaxation times and the initial shear rate or stress in step tests, which improves model predictions in these tests. The best-fit value of

*n*depends on the choice of

*ϕ*(i.e., stress or shear rate) and can be experimentally determined from the fast response in shear start-up tests at large shear rates or stress [33]. In our previous work, Eq. (4) was written in either a stress-controlled (SC) form (

*ϕ*=

*σ*) or a rate-controlled (RC) one (

*ϕ*= $\gamma \u0307$). In the ML-IKH model, with the assumption that thixotropy occurs only upon plastic deformation,

*ϕ*is written as follows:

The above equation states that *λ* does not decrease within the yield surface where $\gamma \u0307p=0$ and $|\sigma eff|\u2264\lambda ky$, and changes in the elastic deformation do not lead to the breakage of internal structures. It is the large and irreversible deformation (i.e., the plastic deformation) that causes thixotropic response. This distinction between reversible elastic strains and irreversible plastic ones, and the dependence of thixotropy only on the latter, are also adopted in the IKH and MDT models.

Note that the RC and SC forms are two similar but different models. They make different constitutive assumptions about how flow conditions influence the thixotropic kinetics. We refer the readers to [55] for a conceptual discussion about the choice of flow parameters, and [56] for a numerical comparison between different choices. A third option based on energy consideration has also been adopted [57,58].

A mechanical analog can be drawn for the proposed model, as shown in Fig. 2. Prior to yielding, $\gamma \u0307p=0$ and this model behaves as a Kelvin-Voigt viscoelastic solid. Upon yielding, the model predicts a Jeffreys-type viscoelastic behavior with a yield stress in unidirectional shearing given by $\lambda ky+khA$ and in the opposite direction of $\u2212\lambda ky+khA$. As *λ* decreases, the yield surface shrinks and the viscoelastic relaxation time ($\lambda \eta thi/G$) decreases. In limit of $\gamma \u0307=+\u221e$ and $\lambda =0$, this model reduces to a Newtonian fluid with a viscosity of $\eta m$.

### B. Tensorial form

In this section, we present the tensorial formulation of the proposed model that can be used in arbitrary 3D flow histories. We write the model in the Eulerian form based on an additive decomposition of the rate of deformation tensor. Compared with the Lagrangian form that is based on the multiplicative decomposition of the deformation gradient tensor, and therefore requires the definition of a fixed reference state with respect to which total strain is defined, the Eulerian form is more suitable for most fluid flow simulations of thixotropic materials. In addition, we assume that the elastic deformation is small while the plastic deformation and rotation might be large. This assumption greatly simplifies the tensorial constitutive equations and is valid for weakly aggregated thixotropic suspensions in which the internal structures can withstand only small elastic strains before breaking.

As in Eq. (11), the total extra stress tensor (**τ**) is divided into contributions from the medium (**τ**_{m}) and from the internal structure (**τ**_{s}), with **τ**_{m} proportional to the rate of deformation tensor (**D**)

Here, **D** is the sum of its elastic and plastic components (**D**_{e} and **D**_{p})

Below, Eq. (25) expresses **τ**_{s} in a differential and frame-invariant form. “▽” denotes the upper-convected time derivative

Its operation on an arbitrary tensor, **S**, is defined in Eq. (26). Here, the overdot denotes the material time derivative, the superscript “T” means transpose, and $\u2207v$ is the velocity gradient tensor

The effective stress tensor (**τ**_{eff}), the counterpart of $\sigma eff$ in the 1D model, is defined as

where $\kappa back$ is an intermediate variable that depends on the back stress tensor ($\tau back=khA$) according to Eq. (28).

Here, **A** is the 3D generalization of the scalar kinematic hardening parameter, *A*. $\kappa back$ is introduced to avoid the aphysical oscillation of **A** in simple shear flows (see details in Sec. VI B). This method was first proposed by Tsakmakis [59] in his finite plasticity model.

Assuming the von Mises yield criterion, the 3D generalization of the plastic flow rule in Eqs. (15)–(18) are written as follows:

where $\sigma \xaf$ is the equivalent shear stress, defined as

The superscript “D” means the deviatoric part of the tensor; e.g., $\tau effD=\tau eff\u2212trace(\tau eff)\delta /3$, where $\delta $ is the unit tensor. The “:” denotes the scalar product (or double dot product) of two tensors.

The tensorial equivalent of Eq. (19) is written as

where again $\theta =\eta thi/G$ is the viscoelastic relaxation time in the fully structured state. The evolution equation for **A** is given by

Here, **A** is symmetric and equals to the zero tensor in the initial undeformed state. This kinetic equation satisfies the second law of thermodynamics and does not bring aphysical oscillation in simple shear flows (see Sec. VI B for further discussion). The “o” over **A** denotes the Jaumann derivative. Its operation on an arbitrary tensor, $S$, is defined as

Here, $W=(\u2207v\u2212(\u2207v)T)/2$ is the vorticity tensor. *d*_{p} is the magnitude of **D**_{p}, defined as

The expression for $\u2009\lambda $ and the kinetic equation for $\lambda i$ are the same as those in the scalar model except that the flow parameter *ϕ* is now written as a function of $\sigma \xaf$ or *d*_{p}, as shown in Eq. (35). *ϕ* is therefore a scalar isotropic function of **τ**_{eff} or **D**_{p} and therefore independent of an arbitrary rotation imposed on the frame.

### C. Summary of the proposed model

Table I is a summary of the proposed model, which has 12 parameters. (In comparison, the thixotropic only ML model has nine parameters, the IKH model has eight, and the MDT model has 13, cf. Table III.) The parameters governing the relationship between stress and shear rate at steady state are $a$, $n$, $kh$, $\u2009ky$, $k1/k3$, $k2/k3$, $\eta thi$, and $\eta m$.

Scalar constitutive equations . | Tensorial constitutive equations . |
---|---|

$\sigma =\eta m\gamma \u0307+\sigma s$ | $\tau =2\eta mD+\tau s$ |

$\lambda \theta \sigma \u0307s+max(0,|\sigma eff|\u2212\lambda ky|\sigma eff|)\sigma eff=\lambda \eta thi\gamma \u0307$ | $\lambda \theta \u2207\tau s+max{0,\sigma \xaf\u2212\lambda ky\sigma \xaf}\tau eff=2\lambda \eta thiD$ |

$\sigma eff=\sigma s\u2212\sigma back$ | $\tau eff=\tau s\u2212\kappa back$ |

$\sigma back=khA$ | $\kappa back=\tau back+2\tau back2/kh,\u2009\tau back=khA$ |

$\gamma \u0307p={0if\u2009|\sigma eff|<\lambda ky|\sigma eff|\u2212\lambda ky\lambda \eta thi\u22c5sign(\sigma eff)if\u2009|\sigma eff|\u2265\lambda ky$ | $Dp={0if\u2009\u2009\u2009\sigma \xaf<\lambda ky(\sigma \xaf\u2212\lambda ky)2\lambda \eta thi\u22c5\tau eff\sigma \xafif\u2009\u2009\u2009\sigma \xaf\u2265\lambda ky$ |

$A\u0307=\gamma \u0307p\u2212qA|\gamma \u0307p|$ | $\omicron A=Dp\u22c5A+A\u22c5Dp+Dp\u2212qdpA$ |

${SC\u2009form:\u2009\varphi =max(0,|\sigma eff|\u2212\lambda ky)RC\u2009form:\u2009\varphi =|\gamma \u0307p|\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009$ | ${SC\u2009form:\u2009\varphi =max(0,\sigma \xaf\u2212\lambda ky)RC\u2009form:\u2009\varphi =2dp$ |

$\lambda =\u2211i=1NCi\lambda i$ | |

$d\lambda idt=Di[\u2212k1\varphi a\lambda in+k2\varphi 0.5(1\u2212\lambda i)+k3(1\u2212\lambda i)]$ | |

$exp\u2009(\u2212t\beta )\u2248\u2211i=1NCi(\beta )\u22c5e\u2212Dit$ | |

Model parameters (12 in total): $governing\u2009steady\u2009state:\u2009a,n,ky,khq,k1k3,k2k3,\eta m,\eta thicontrolling\u2009transient\u2009response:\u2009k3,\beta ,q,\u2009G$ |

Scalar constitutive equations . | Tensorial constitutive equations . |
---|---|

$\sigma =\eta m\gamma \u0307+\sigma s$ | $\tau =2\eta mD+\tau s$ |

$\lambda \theta \sigma \u0307s+max(0,|\sigma eff|\u2212\lambda ky|\sigma eff|)\sigma eff=\lambda \eta thi\gamma \u0307$ | $\lambda \theta \u2207\tau s+max{0,\sigma \xaf\u2212\lambda ky\sigma \xaf}\tau eff=2\lambda \eta thiD$ |

$\sigma eff=\sigma s\u2212\sigma back$ | $\tau eff=\tau s\u2212\kappa back$ |

$\sigma back=khA$ | $\kappa back=\tau back+2\tau back2/kh,\u2009\tau back=khA$ |

$\gamma \u0307p={0if\u2009|\sigma eff|<\lambda ky|\sigma eff|\u2212\lambda ky\lambda \eta thi\u22c5sign(\sigma eff)if\u2009|\sigma eff|\u2265\lambda ky$ | $Dp={0if\u2009\u2009\u2009\sigma \xaf<\lambda ky(\sigma \xaf\u2212\lambda ky)2\lambda \eta thi\u22c5\tau eff\sigma \xafif\u2009\u2009\u2009\sigma \xaf\u2265\lambda ky$ |

$A\u0307=\gamma \u0307p\u2212qA|\gamma \u0307p|$ | $\omicron A=Dp\u22c5A+A\u22c5Dp+Dp\u2212qdpA$ |

${SC\u2009form:\u2009\varphi =max(0,|\sigma eff|\u2212\lambda ky)RC\u2009form:\u2009\varphi =|\gamma \u0307p|\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009\u2009$ | ${SC\u2009form:\u2009\varphi =max(0,\sigma \xaf\u2212\lambda ky)RC\u2009form:\u2009\varphi =2dp$ |

$\lambda =\u2211i=1NCi\lambda i$ | |

$d\lambda idt=Di[\u2212k1\varphi a\lambda in+k2\varphi 0.5(1\u2212\lambda i)+k3(1\u2212\lambda i)]$ | |

$exp\u2009(\u2212t\beta )\u2248\u2211i=1NCi(\beta )\u22c5e\u2212Dit$ | |

Model parameters (12 in total): $governing\u2009steady\u2009state:\u2009a,n,ky,khq,k1k3,k2k3,\eta m,\eta thicontrolling\u2009transient\u2009response:\u2009k3,\beta ,q,\u2009G$ |

The scalar model predicts that the values of $A$, $\lambda $, and $\sigma $ at steady state are as follows:

The function *f* in Eq. (37) represents the solution of Eq. (4) at steady state where $d\lambda i/dt=0$ for all $\lambda i$.

The parameters governing the transient response are $k3$, $\beta $, $q$, and $G$. $k3$ directly influences the time scale of thixotropic relaxation. $\beta $ determines the extent of deviation from single exponential relaxation. For $\beta =1$, the stretched-exponential reduces to a single exponential and only a single $\lambda i$ is needed, which reduces the model to a single-lambda model. $q$ determines the plastic strain at which $A$ reaches steady state. $\eta thi/G$ governs the time scale of the viscoelastic relaxation.

The model reduces to several different models in certain limits, as described in Table II.

Model parameters . | $\theta $ . | $kh$ . | $ky$ . | $k1,k2,\u2009and\u2009k3$^{a}
. | $\beta \u22121$ . | $\eta m$ . |
---|---|---|---|---|---|---|

Oldroyd-B model | = 0 | = 0 | = 0 | — | ||

Saramito's model [51] | = 0 | = 0 | — | |||

Upper-Convected Maxwell model | = 0 | = 0 | = 0 | — | = 0 | |

ML model in [33] | = 0 | = 0 | ||||

IKH model^{b} | = 0 | = 0 | ||||

Monolambda ideal thixotropic models | = 0 | = 0 | = 0 | = 0 | ||

Bingham model | = 0 | = 0 | = 0 | — | ||

Newtonian fluid | = 0 | = 0 | = 0 | = 0 | — | |

Implication of zero values | No elasticity | No back stress | No yield stress | No thixotropic model component | Single $\lambda $ |

Model parameters . | $\theta $ . | $kh$ . | $ky$ . | $k1,k2,\u2009and\u2009k3$^{a}
. | $\beta \u22121$ . | $\eta m$ . |
---|---|---|---|---|---|---|

Oldroyd-B model | = 0 | = 0 | = 0 | — | ||

Saramito's model [51] | = 0 | = 0 | — | |||

Upper-Convected Maxwell model | = 0 | = 0 | = 0 | — | = 0 | |

ML model in [33] | = 0 | = 0 | ||||

IKH model^{b} | = 0 | = 0 | ||||

Monolambda ideal thixotropic models | = 0 | = 0 | = 0 | = 0 | ||

Bingham model | = 0 | = 0 | = 0 | — | ||

Newtonian fluid | = 0 | = 0 | = 0 | = 0 | — | |

Implication of zero values | No elasticity | No back stress | No yield stress | No thixotropic model component | Single $\lambda $ |

^{a}

In addition, $\lambda $ is set to be unity.

^{b}

With additional requirements: $k2=0$, $a=1$, $n=1$, $\varphi =|\gamma \u0307p|$, replacing $\lambda \eta thi$ with $\eta thi$, and using a different evolution equation for *A* or **A** (see detailed comparison in Sec. VI B).

### D. Qualitative comparison with the ML, IKH, and MDT models

We first qualitatively compare the ML-IKH model with the ML, IKH, and MDT models. Table III is a summary. As reviewed in Sec. II, TEVP materials behave like viscoelastic solids before yielding and viscoelastic fluids after yielding. Moreover, their response approaches ideal thixotropy at high shear rates (where the required shear rate level may differ in different systems. For the fumed silica suspensions reported in [24,33,40], $\gamma \u0307>5\u2009s\u22121$). The thixotropic relaxation has multiple time scales and follows a stretched exponential in step tests. In general step-shear-rate tests, the stress exhibits an instant jump, followed by a fast elastic relaxation with an initial extremum and a subsequent slow thixotropic relaxation. The ML, IKH, and MDT models predict only partially these rheological behaviors. All four models fail to predict the delayed yielding behavior that some TEVP materials exhibit. We refer the readers to [43,60–62] for models that address the delayed yielding phenomenon.

. | ML . | IKH . | MDT . | ML-IKH . |
---|---|---|---|---|

Number of parameters | 9 | 8 | 13 | 12 |

Preyield response | Rigid solid | Linear elastic solid | Linear elastic solid | Linear viscoelastic solid |

Response in fully unstructured state | Newtonian fluid | Maxwell fluid | Newtonian fluid | Newtonian fluid |

Kinematic hardening | None | Modified Armstrong-Frederick | None | Armstrong-Frederick |

Prediction of delayed yielding | No | No | No | No |

Predictions in step-rate tests: (1) Instantaneous stress jump | Yes | No | Yes | Yes |

(2) Fast elastic relaxation with an initial extremum | No | Yes | No | Yes |

(3) Stretched-exponential thixotropic relaxation | Yes | No | No | Yes |

Summary of the IKH model $\gamma =\gamma e+\gamma p,\u2009\u2009\u2009\sigma =G\gamma e,\u2009\u2009\u2009A\u0307=\gamma \u0307p\u2212(q|A|)msign(A)|\gamma \u0307p|,\u2009\u2009\u2009\lambda \u0307=k1(1\u2212\lambda )\u2212k3\lambda |\gamma \u0307p|\gamma \u0307p=max{0,sign(|\sigma \u2212khA|)(|\sigma \u2212khA|\u2212\lambda ky)/\mu p}Model\u2009parameters:\u2009G,q,m,k1,k3,kh,ky,\mu p$ | ||||

Summary of the MDT model $\gamma =\gamma e+\gamma p,\u2009\u2009\u2009\gamma \u0307e=\gamma \u0307p\u2212(\gamma e|\gamma \u0307p|)/\gamma max,\u2009\u2009\u2009\gamma max=min(\gamma CO/\lambda m,1),\u2009\u2009\u2009\gamma CO=\sigma y0/G0G\u0307f=\u2212kG(Gf\u2212\lambda G0),\u2009\u2009\u2009\lambda \u0307=kBrown[\u2212\lambda |t\u0302r1\gamma \u0307p|a+(1\u2212\lambda )(1+|t\u0302r2\gamma \u0307p|d)],\u2009\u2009\u2009\sigma =Gf\gamma e+(\lambda KST|\gamma \u0307p|n2+K\u221e|\gamma \u0307p|n1)sign(\gamma \u0307p)Model\u2009parameters:\u2009m,G0,\sigma y0,kG,kBrown,t\u0302r1,t\u0302r2,a,d,n1,n2,KST,K\u221e$ |

. | ML . | IKH . | MDT . | ML-IKH . |
---|---|---|---|---|

Number of parameters | 9 | 8 | 13 | 12 |

Preyield response | Rigid solid | Linear elastic solid | Linear elastic solid | Linear viscoelastic solid |

Response in fully unstructured state | Newtonian fluid | Maxwell fluid | Newtonian fluid | Newtonian fluid |

Kinematic hardening | None | Modified Armstrong-Frederick | None | Armstrong-Frederick |

Prediction of delayed yielding | No | No | No | No |

Predictions in step-rate tests: (1) Instantaneous stress jump | Yes | No | Yes | Yes |

(2) Fast elastic relaxation with an initial extremum | No | Yes | No | Yes |

(3) Stretched-exponential thixotropic relaxation | Yes | No | No | Yes |

Summary of the IKH model $\gamma =\gamma e+\gamma p,\u2009\u2009\u2009\sigma =G\gamma e,\u2009\u2009\u2009A\u0307=\gamma \u0307p\u2212(q|A|)msign(A)|\gamma \u0307p|,\u2009\u2009\u2009\lambda \u0307=k1(1\u2212\lambda )\u2212k3\lambda |\gamma \u0307p|\gamma \u0307p=max{0,sign(|\sigma \u2212khA|)(|\sigma \u2212khA|\u2212\lambda ky)/\mu p}Model\u2009parameters:\u2009G,q,m,k1,k3,kh,ky,\mu p$ | ||||

Summary of the MDT model $\gamma =\gamma e+\gamma p,\u2009\u2009\u2009\gamma \u0307e=\gamma \u0307p\u2212(\gamma e|\gamma \u0307p|)/\gamma max,\u2009\u2009\u2009\gamma max=min(\gamma CO/\lambda m,1),\u2009\u2009\u2009\gamma CO=\sigma y0/G0G\u0307f=\u2212kG(Gf\u2212\lambda G0),\u2009\u2009\u2009\lambda \u0307=kBrown[\u2212\lambda |t\u0302r1\gamma \u0307p|a+(1\u2212\lambda )(1+|t\u0302r2\gamma \u0307p|d)],\u2009\u2009\u2009\sigma =Gf\gamma e+(\lambda KST|\gamma \u0307p|n2+K\u221e|\gamma \u0307p|n1)sign(\gamma \u0307p)Model\u2009parameters:\u2009m,G0,\sigma y0,kG,kBrown,t\u0302r1,t\u0302r2,a,d,n1,n2,KST,K\u221e$ |

### E. Physical interpretation of structure parameters

In the ML-IKH model, {*λ _{i}*},

*A*, and

*γ*

_{e}account phenomenologically for, respectively, the instantaneous degree of structures (in the isotropic sense), characteristic orientation, and elastic deformation of the internal structures. In single-lambda thixotropic models, the structure parameter $\lambda $ is often interpreted as the instantaneous number of links (or other structural units) normalized by the maximum number in a fully structured state [25,63]. As a first approximation, it is assumed that all links have the same strength and break or form at the same rate. Most thixotropic models assume a first-order mechanism of the breakage, which results in a linear kinetic equation for

*λ*. The use of a set of {

*λ*} is a generalization of such description. We envision {

_{i}*λ*} as the relative numbers of internal links of different types specified by their relative rates of breakage or formation (the

_{i}*D*values). The yield stress and viscosity are related to the linear combination of

_{i}*λ*with the weights following a stretched exponential distribution.

_{i}As discussed in [15,64], *A* quantifies the level of structural anisotropy, or characteristic orientation of internal structures, which gives rise to the back stress. *A = *0 implies randomly orientated microstructures and *A = *1/*q* gives the steady-state anisotropic orientation in simple shear flows. In 3D deformations, the scalar parameter *A* is replaced by the tensor **A**. (Since our model accounts for multiple types of structures through the multiple *λ _{i}* values, it might make sense to include multiple

*A*values, one corresponding to each type of structure, and possibly each having a different value of

*q*. However, for simplicity, we forego this refinement here.) Shear-induced anisotropic microstructures of thixotropic fluids have been experimentally observed in Laponite suspensions [65], fumed silica suspensions [66], waxy crude oils [15], etc. A review is given by Vermant and Solomon [2].

The internal structures withstand a certain amount of elastic deformation before breaking, and the relaxation of this elastic deformation results in the viscoelastic response of thixotropic materials. *γ*_{e} can be interpreted as the characteristic elastic strain of internal structures. One 3D generalization of *γ*_{e} is the tensor strain $(Be\u2212\delta )$ where $Be$ is the elastic Finger tensor, defined as $Be=Fe\u22c5FeT$. Here, $Fe$ is the elastic component of the deformation gradient (**F**) according to the Kroner Decomposition [67], $F=Fe\u22c5Fp$, where $Fp$ is the plastic deformation gradient. A detailed discussion about this decomposition is provided in the Appendix.

## IV. MODEL PARAMETERIZATION

The parameterization of the proposed model requires the following experimental data: The steady state flow curve, a flow reversal test with small shear rate, and several step-shear-rate tests with relatively small initial shear rates. For example, in this study, we use the experimental data in Fig. 3 and the flow reversal test with $\gamma \u03070=0.01\u2009s\u22121$ in Fig. 6(a) to parameterize the model.

As discussed before, the model predicts a viscosity of *η*_{m} in the fully unstructured state and a yield stress of (*k*_{h}/*q + k*_{y}) in the totally structured state. Therefore, *η*_{m} can be directly estimated from the high-shear-rate limit of the steady state flow curve and (*k*_{h}/*q + k*_{y}) from the zero-shear-rate limit.

In a flow reversal test where the shearing direction instantaneously changes while the shear rate remains the same, the model predicts a two-step relaxation of shear stress. This two-step relaxation is caused by the evolution of *γ*_{e}, the viscoelastic relaxation, and that of *A*, the kinematic hardening relaxation. In the case where the shear rate is small and viscous stress is insignificant, the amount of stress that undergoes the kinematic hardening relaxation, as predicted by this model, equals to *k*_{h}/*q*, which can be estimated experimentally. Next, the value of *q*, governing the plastic strain at which *A* reaches steady state in a flow reversal test, can be also determined in the same flow reversal test through best fit.

The remaining parameters are divided into two groups—one governing the steady state response, including *a*, *n*, $k1/k3$, $k2/k3$ and *η*_{thi}; the other one governing the transient response, including *k*_{3}, *β*, and *G*. Parameters in the first group are first determined from the best fit of the steady state flow curve and those in the second group from the best fit of several step-shear-rate tests [in this study, we use three step tests from 0.1 s^{–1}, as shown in Fig. 3(b)]. Iteration may be needed to adjust the boundaries of the parameters in the first group, particularly $k1/k3$, $k2/k3,$ and *η*_{thi}. Alternatively, the power-law exponent *n* in the shear-induced structure break-up term can be determined from the fast response in shear startup tests at large shear rates or stresses. In such tests, the viscous stress dominates and the model predicts an apparent viscosity of $(\lambda \eta thi+\eta m)$. Therefore, the transients of *λ* can be estimated from the measured viscosity. In the very short period after the stepwise change of shear rate or stress, if $n=1$, the model predicts a linear relationship between $ln\lambda $ and time; while if $n\u22601$, a linear relationship is predicted between $\lambda 1\u2212n$ and time. The value of $n$ may be different for the SC and RC forms. More details about the estimation of $n$ are given in our previous work [33].

The ${Ci}$ and ${Di}$ values as well as *N*, the number of substructure parameters, are determined from the value of *β*. As shown in Eq. (21), they are calculated to best approximate the stretched exponential using a limited number of mono-exponentials. More details including a matlab subroutine to perform this calculation are provided in the supplementary material [39].

## V. MODEL EVALUATION

In this section, we test the proposed model, in both SC and RC forms, against experimental results with shear histories including steady state, step shear rate, step stress, intermittent step shear rate, flow reversal, and LAOS. We also investigate the predictions of the ML, IKH, and MDT models for comparison. There are other models in the literature that fall into the category of TEVP models. One important example is the “unified-approach” (UA) model of de Souza Mendes and Thompson [43]. It has been shown that the UA model qualitatively captures the TEVP rheology of the fumed silica suspension but lacks quantitative agreement [24]. For brevity, we do not include it in this work.

We discuss the critical model components and corresponding improvements in predicting each type of tests. The LAOS data are adopted from the work of Armstrong *et al*. [24] and the others from our previous work [33]. All rheological tests use cone and plate geometries. Both data sets are for fumed silica suspensions (a hydrophobic fumed silica, Aerosil^{®} R972, is dispersed in a blend of paraffin oil and low-molecular-weight polymer). The two fumed silica suspensions were synthesized according to the same formula reported in [23] except that they contain different polymer ingredients. For more details about the sample preparation and experimental protocols, we refer the readers to [24,33].

Note that in this section the original IKH model (summarized in Table III) is modified by replacing the constant plastic viscosity, $\mu p$, with $\lambda \eta thi+\eta m$. This modification is necessary because the original IKH model predicts that thixotropic response vanishes when the viscous stress dominates, which does not agree with the rheology of the fumed silica suspension [24,33,37,40] (e.g., the original IKH model would predict a nearly constant viscosity of $\mu p$ in the tests of Fig. 4, contradicting the experimental results). This modification introduces one additional parameter. In the following discussion, we refer to this slightly modified model as IKH-V model. The predictions of the original IKH model are provided in the supplementary material [39].

### A. Steady state flow curve

Figure 3(a) shows the steady state flow curve of the fumed silica suspension reported in [33] and the fitting results of the ML-IKH model in both SC and RC forms, the IKH-V, MDT, and ML (in its SC form) models. The steady state flow curve can be fitted accurately by all models.

Table IV lists the parameter values that generate the model calculations. The MDT model is parameterized from the best fit of the experimental data in Fig. 3. Besides these data, the parameterization of the IKH-V and ML-IKH models also requires the experimental results of the flow reversal test with a shear rate of 0.01 s^{–1} in Fig. 6. More details about the parameterization of the MDT and IKH models are provided in the supplementary material [39].

Parameters . | SC form . | RC form . | ML(SC) . | IKH-V model . | MDT model . |
---|---|---|---|---|---|

$ky$ | $0.434\u2009Pa$ | $0.458\u2009Pa$ | 0.68 Pa | $0.51\u2009Pa$ | $\sigma y0=0.8Pa$, $K\u221e=0.41Pa\u22c5s$, $n1=1$, $G0=100\u2009Pa$, $m=0.52$, $n2=1$, $KST=1.76\u2009Pa\u2009s$, $a=0.7$, $d=0.3$, $kBrown=8.7\xd710\u22122\u2009s\u22121$, $t\u0302r1=3.42\u2009s$, $t\u0302r1=0.18\u2009s$, $kG=8.2\xd710\u22122\u2009s\u22121$ |

$kh$ | $0.195\u2009Pa$ | $0.205\u2009Pa$ | 0 | $0.22\u2009Pa$ | |

$q$ | 1 | 1 | 0 | $1$ | |

$G$ | $60\u2009Pa$ | $70\u2009Pa$ | +$\u221e$ | $100\u2009Pa$ | |

$\eta thi$ | $2.325\u2009Pa\u2009\u22c5s$ | $1.840\u2009Pa\u2009\u22c5s$. | 1.74 $Pa\u2009\u22c5s$ | $1.29\u2009Pa\u2009\u22c5s$ | |

$\eta m$ | $0.430\u2009Pa\u2009\u22c5s$ | $0.425\u2009Pa\u2009\u22c5s$ | 0.45 $Pa\u2009\u22c5s$ | $0.56\u2009Pa\u2009\u22c5s$ | |

$\beta $ | 0.65 | 0.65 | 0.65 | ||

$a$ | 2.537 | 1.758 | 1.9 | ||

$n$ | 2 | 3 | 2 | ||

$k1$ | $0.042$$Pa\u2212as\u22121$ | $0.105$ $sa\u22121$ | 0.043 $Pa\u2212as\u22121$ | $6.5\xd710\u22122$ | |

$k2$ | $3\xd710\u22124\u2009Pa\u22120.5s\u22121$ | $3\xd710\u22124\u2009s\u22120.5$ | 0.0104 $Pa\u22120.5s\u22121$ | $7.4\xd710\u22122\u2009s\u22121$ | |

$k3$ | $6\xd710\u22123\u2009s\u22121$ | $1.6\xd710\u22122\u2009s\u22121$ | 0.01 $s\u22121$ | $m=1$ |

Parameters . | SC form . | RC form . | ML(SC) . | IKH-V model . | MDT model . |
---|---|---|---|---|---|

$ky$ | $0.434\u2009Pa$ | $0.458\u2009Pa$ | 0.68 Pa | $0.51\u2009Pa$ | $\sigma y0=0.8Pa$, $K\u221e=0.41Pa\u22c5s$, $n1=1$, $G0=100\u2009Pa$, $m=0.52$, $n2=1$, $KST=1.76\u2009Pa\u2009s$, $a=0.7$, $d=0.3$, $kBrown=8.7\xd710\u22122\u2009s\u22121$, $t\u0302r1=3.42\u2009s$, $t\u0302r1=0.18\u2009s$, $kG=8.2\xd710\u22122\u2009s\u22121$ |

$kh$ | $0.195\u2009Pa$ | $0.205\u2009Pa$ | 0 | $0.22\u2009Pa$ | |

$q$ | 1 | 1 | 0 | $1$ | |

$G$ | $60\u2009Pa$ | $70\u2009Pa$ | +$\u221e$ | $100\u2009Pa$ | |

$\eta thi$ | $2.325\u2009Pa\u2009\u22c5s$ | $1.840\u2009Pa\u2009\u22c5s$. | 1.74 $Pa\u2009\u22c5s$ | $1.29\u2009Pa\u2009\u22c5s$ | |

$\eta m$ | $0.430\u2009Pa\u2009\u22c5s$ | $0.425\u2009Pa\u2009\u22c5s$ | 0.45 $Pa\u2009\u22c5s$ | $0.56\u2009Pa\u2009\u22c5s$ | |

$\beta $ | 0.65 | 0.65 | 0.65 | ||

$a$ | 2.537 | 1.758 | 1.9 | ||

$n$ | 2 | 3 | 2 | ||

$k1$ | $0.042$$Pa\u2212as\u22121$ | $0.105$ $sa\u22121$ | 0.043 $Pa\u2212as\u22121$ | $6.5\xd710\u22122$ | |

$k2$ | $3\xd710\u22124\u2009Pa\u22120.5s\u22121$ | $3\xd710\u22124\u2009s\u22120.5$ | 0.0104 $Pa\u22120.5s\u22121$ | $7.4\xd710\u22122\u2009s\u22121$ | |

$k3$ | $6\xd710\u22123\u2009s\u22121$ | $1.6\xd710\u22122\u2009s\u22121$ | 0.01 $s\u22121$ | $m=1$ |

### B. Step tests

Step tests refer to the experiments where the shear rate or stress undergoes a stepwise change after a period of steady shear. The transient response in step tests results mainly from thixotropy and viscoelasticity.

Figure 3(b) shows the experimental results and model fits in three step-shear-rate tests from 0.1 s^{–1} to different final shear rates. The stress exhibits a fast elastic relaxation with an initial maximum followed by a gradual thixotropic relaxation. The thixotropic response follows a stretched-exponential relaxation. The ML(SC) and ML-IKH models successfully capture this behavior through the ML approach. The IKH-V model and MDT model, however, always show a single exponential thixotropic relaxation, which results in the sharper transition of stress transients than the experimental results. The IKH-V and ML-IKH models fit the initial increase and maximum of stress while the ML(SC) and MDT models do not. The SC form of the ML-IKH model yields slightly better fitting results than does the RC form.

In step tests with relatively large initial shear rates or stress, the fast elastic response is insignificant, and the test material behaves like an ideal thixotropic fluid. While the stress in step-shear-rate tests shows an instantaneous jump, the measured viscosity changes continuously. In this region of ideal thixotropy, the thixotropic-only ML model and the ML-IKH model have a similar level of agreement with experimental results, as shown in Fig. 4. (The SC and RC forms of the ML-IKH model have very similar predictions, as shown in Fig. 11. For brevity, in Fig. 4, we show the predictions of only the SC form.)

Figure 4(a) shows the model predictions in a step-shear-rate test from 3 to 9 s^{–1}. Compared with their results in Fig. 3(b), the IKH-V model and MDT model show better predictions in this test. The MDT model agrees with experimental results well although it still shows a monoexponential stress relaxation. The IKH-V model can only qualitatively capture the thixotropic response. The ML and ML-IKH models slightly overpredict the stress, but it agrees with the experimental results better than do the other two models, due to the stretched-exponential relaxation in it.

The results of two step-stress tests (from 3 and 9 Pa to 7 Pa) are given in Fig. 4(b). Instrument inertia causes the steep response for times less than 0.1 s, and is accounted for theoretically as described in our previous work [33]. Both the IKH-V and MDT models show significant deviations from the experimental data while the ML and ML-IKH models give more accurate predictions.

Next, we examine the models in more step tests with different test conditions, and present these results in Appendix A, Fig. 11. This figure, along with Figs. 3(b) and 4 show that the ML-IKH model can quantitatively predict a wide range of step tests while the MDT model and IKH-V model fail to predict as accurately over this range. In the ideal thixotropy region, multiple thixotropic structure parameters (corresponding model parameter: *β*) and their nonlinear kinetic equations (*a*, *n*) are two key factors needed to predict accurately these step tests. For step tests with small initial shear rates, viscoelasticity (*G*) is important in capturing the fast elastic response. The AF equation (*k*_{h}, *q*, *m*) does not contribute to the transient response because the shearing direction remains constant.

### C. Intermittent step tests

Intermittent step tests have the flow histories of two consecutive step tests where the duration of the first one is varied. Figure 5(a) shows the shear rate histories of a group of intermittent step tests. Corresponding experimental data and model predictions are shown in Figs. 5(b)–5(f). The stress transients for time >0 are sensitive to the value of Δ*t* up to about 10 min. The stress evolves monotonically when Δ*t* is small or large while for some intermediate values of Δ*t*, nonmonotonic *thixotropic* evolution (stress transients for time >0.1 s) occurs—stress first decreases, then builds up, and eventually reaches steady state.

Such nonmonotonic thixotropic relaxations should be distinguished from the elastic response, which has a much shorter time scale (<0.1 s). The ML Approach successfully captures this behavior. The rate of change of each *λ _{i}* is proportional to the prefactor

*D*. In the test of Fig. 5(a), initially all

_{i}*λ*have the same steady state values at 100 s

_{i}^{–1}. After the shear rate steps down to 0.5 s

^{–1}, the {

*λ*} build up at different rates. For intermediate Δ

_{i}*t*values, the

*λ*values with large

_{i}*D*build up enough that they then decrease as the shear rate steps to 5 s

_{i}^{–1}. The

*λ*with small

_{i}*D*, however, have not rebuilt much and so continue to increase after the stepwise change of shear rate. Consequently, in the period of

_{i}*t*> 0, the stress first decreases and then increases.

The IKH-V and MDT models incorporate a single structure parameter for thixotropic evolution and therefore show poor agreement with experimental results, as shown in Figs. 5(b) and 5(c). Figure 5(d) shows that the ML-IKH model captures both the nonmonotonic thixotropic evolution and the dependence of the stress transients on Δ*t*. (For brevity, we show only the predictions of SC form. The RC form has similar level of agreement). In Fig. 5(e), we show the predictions of the ML-IKH model with the same parameters values except that *β* is set to unity, so that *N* = *D _{i}* =

*C*1, and the ML-IKH model reduces to a single-lambda model. Its predictions in this limit, like those of the IKH-V and MDT models, show a monotonic

_{i}=*thixotropic*relaxation when the shear rate is constant. And the overall predictions are worse than those in Fig. 5(c). Note that the initial elastic response is only slightly influenced by the value of

*β*. The elastic response might also lead to nonmonotonic stress transients but the elastic response has a much smaller time scale than the thixotropic relaxation. Figure 5(f) shows the predictions of the thixotropic-only ML model. The ML and ML-IKH models have a similar level of agreement with experimental results, which indicates that the ML approach is the dominating model component in predicting the test of Fig. 5(a). For shear histories that involve further smaller shear rates than 0.5 s

^{–1}(see Fig. S5 for an example), the differences between ML and ML-IKH model predictions increase. The latter gives better predictions because viscoelastic evolution is significant when low shear rates are involved and the ML model lacks the component accounting for viscoelasticity.

In summary, the ML Approach (corresponding model parameter, *β*), and the viscoelastic model component (*G*) if low shear rates are involved, are critical to qualitatively predict intermittent step tests. The ML approach successfully predicts the nonmonotonic *thixotropic* evolution. The effects of nonlinear thixotropic evolution equation (*a*, *n*) cannot be clearly seen here. The AF equation (*k*_{h}, *q*, *m*) does not contribute to the transient response because the shearing direction remains the same.

### D. Flow reversal tests

In a flow reversal test, after an attainment of a steady state under shearing, an instantaneous reversal of the shear direction is imposed, with the magnitude of the shear rate held fixed.

For a Newtonian fluid in a flow reversal test, the sign of the shear stress instantaneously changes while its magnitude remains constant. A Maxwell fluid, on the other hand, predicts that the stress after the reversal continuously evolves to its opposite value following an exponential function with a constant characteristic time. For some colloidal suspensions that form anisotropic structures under shear, the stress following flow reversal exhibits a gradual relaxation that scales with strain [68,69]. TEVP fluids, e.g., the thixotropic fumed silica suspension in this study, however, show different rheological behaviors, as shown in Fig. 6.

In Fig. 6, we compare the experimental results with model fits (for $\gamma \u03070=0.01\u2009s\u22121$) or independent predictions (for all $\gamma \u03070>0.01\u2009s\u22121$) in several flow reversal tests, each identified by a specific color of symbols and lines. Figure 6(a) shows experimental data plotted as reduced stress versus strain in flow reversal tests with shear rates ranging from 0.01 to 10 s^{–1}. The data were collected using the strain-controlled ARES rheometer. The test fluid and procedures of rheological treatment are as previously reported [33]. Similar results have been reported in [70]. In flow reversal tests with shear rates of 0.1 s^{–1} or lower, the reduced stress can be scaled with the strain to yield a master curve. The stress transients follow a two-step relaxation that takes approximately three strain units to reach the plateau at steady state. As the shear rate increases above 0.1 s^{–1}, the initial steep response is no longer detectable and the reduced stress transients deviate from the master curve.

Figure 6(b) shows the model fits and predictions corresponding to the experimental results of Fig. 6(a). The flow reversal test with $\gamma \u03070=0.01\u2009s\u22121$ is used to fix the values of *k*_{h}, *q*, and *m*. The optimal values of *q* and *m* are both unity in the IKH-V and ML-IKH models. (The RC and SC forms give similar predictions. For brevity, we show only the predictions of the SC form.) The two models yield similar results in all the flow reversal tests of Fig. 6(a), both predicting the two-step relaxation in good agreement with experimental data. In Fig. 6(b), for clarity, we show only two typical predictions of the IKH-V model.

Since the absolute value of the shear rate remains constant in a flow reversal test, the predicted thixotropic effect is insignificant. The transient response results mainly from the viscoelastic relaxation and the plastic evolution (for example, kinematic hardening effects). In the IKH-V and ML-IKH models, the evolution of *γ*_{e}, the viscoelastic relaxation, captures the initial steep response; and the transients in *A*, which describe the kinematic hardening effects, allow both models to predict well the subsequent gradual stress response that scales with strain. The two models predict the rise of the initial plateau as shear rate increases. But they overpredict the initial plateau and underpredict the amount of strain it takes for the stress to reach steady state. The discrepancies are significant for large shear rates, possibly because the viscosity is also anisotropic and evolves in reversing flows with a smaller time (or strain) scale than that of the back-stress evolution. The MDT model captures only the elastic effects and therefore has poor predictions in flow reversal tests. The ML model lacks components of both viscoelasticity and kinematic hardening. It, therefore, predicts no transient response in flow reversal tests.

### E. LAOS

In this section, we test the ML, ML-IKH, IKH-V, and MDT models in LAOS tests where both the shearing direction and shear rate are periodically changing. The strain history of a LAOS test is a sinusoidal wave with a frequency of *ω* and amplitude of $\gamma 0$.

The experimental data are originally reported in [24]. The test fluid is a fumed silica suspension and the rheological tests were performed on the ARES-G2 rheometer. To parameterize the ML, ML-IKH, and IKH-V models, we used the following experimental data from [24]: The steady state flow curve, the step-shear-rate tests from 5 and 0.1 s^{–1}, and the flow reversal test from −1 s^{−1}. The predictions of the MDT model are calculated using the original parameter values given in [24]. Table V lists all parameter values. For this set of experimental data, the RC forms of the ML and ML-IKH models have better agreement with experimental data than the SC forms. For brevity, in Figs. 7 and 8, we show predictions of only the RC forms of the ML and ML-IKH models.

Parameters . | ML-IKH (RC form) . | ML (RC form) . | IKH-V . | MDT . |
---|---|---|---|---|

$ky$ | $7.66\u2009Pa$ | 10.96 Pa | $ky=7.88\u2009Pakh=3.36q=1m=1G=700\u2009Pa\eta thi=16.32\u2009Pa\u22c5s\eta m=1.32\u2009Pa\u22c5sk1=1.18k2=0.36\u2009s\u22121$ | $\sigma y0=11\u2009PaK\u221e=1.17\u2009Pa\u22c5sn1=1G0=450\u2009Pam=\u22121.5n2=0.81KST=11.21\u2009Pa\u22c5sa=1.53d=0.63kBrown=0.28\u2009s\u22121t\u0302r1=0.77\u2009st\u0302r1=2.0\u2009skG=0.09\u2009s\u22121$ Values are adopted from [24]. |

$kh$ | $3.30\u2009Pa$ | 0 | ||

$q$ | 1 | 0 | ||

$G$ | $700\u2009Pa$ | $+\u221e$ | ||

$\eta thi$ | $20.0\u2009Pa\u22c5s$ | 16.4 $Pa\u22c5s$ | ||

$\eta m$ | $1.22\u2009Pa\u22c5s$ | 1.33 $Pa\u22c5s$ | ||

$\beta $ | 0. 5 | 0.5 | ||

$a$ | 1.29 | 1.45 | ||

$n$ | 1 | 1 | ||

$k1$ | $0.496\u2009sa\u22121$ | 0.24 $sa\u22121$ | ||

$k2$ | $0.178\u2009s\u22120.5$ | 0.2 $s\u22120.5$ | ||

$k3$ | $0.2\u2009s\u22121$ | 0.1 $s\u22121$ |

Parameters . | ML-IKH (RC form) . | ML (RC form) . | IKH-V . | MDT . |
---|---|---|---|---|

$ky$ | $7.66\u2009Pa$ | 10.96 Pa | $ky=7.88\u2009Pakh=3.36q=1m=1G=700\u2009Pa\eta thi=16.32\u2009Pa\u22c5s\eta m=1.32\u2009Pa\u22c5sk1=1.18k2=0.36\u2009s\u22121$ | $\sigma y0=11\u2009PaK\u221e=1.17\u2009Pa\u22c5sn1=1G0=450\u2009Pam=\u22121.5n2=0.81KST=11.21\u2009Pa\u22c5sa=1.53d=0.63kBrown=0.28\u2009s\u22121t\u0302r1=0.77\u2009st\u0302r1=2.0\u2009skG=0.09\u2009s\u22121$ Values are adopted from [24]. |

$kh$ | $3.30\u2009Pa$ | 0 | ||

$q$ | 1 | 0 | ||

$G$ | $700\u2009Pa$ | $+\u221e$ | ||

$\eta thi$ | $20.0\u2009Pa\u22c5s$ | 16.4 $Pa\u22c5s$ | ||

$\eta m$ | $1.22\u2009Pa\u22c5s$ | 1.33 $Pa\u22c5s$ | ||

$\beta $ | 0. 5 | 0.5 | ||

$a$ | 1.29 | 1.45 | ||

$n$ | 1 | 1 | ||

$k1$ | $0.496\u2009sa\u22121$ | 0.24 $sa\u22121$ | ||

$k2$ | $0.178\u2009s\u22120.5$ | 0.2 $s\u22120.5$ | ||

$k3$ | $0.2\u2009s\u22121$ | 0.1 $s\u22121$ |

The test fluid shows complex rheological behaviors in LAOS tests due to the coexistence of thixotropic, viscoelastic, and plastic response. Their prominence vary depending on the values of $\omega $ and $\gamma 0$. Figure 7 shows the results of a LAOS test with $\omega =0.1\u2009rad/s$ and $\gamma 0=10$. This test has a Deborah number of $De=\eta thi\omega /G\u22483\xd710\u22123$, at which the elastic effects are insignificant. We chose this test condition for detailed discussion because it is suitable to examine the model components accounting for thixotropy and kinematic hardening. We compare the models in more test conditions in Fig. 8.

Figures 7(a) and 7(b) plot, respectively, the elastic (reduced stress vs strain) and viscous (reduced stress vs shear rate) projection of the Lissajous–Bowditch curve. Experimental data and model predictions are both normalized into the range −1−1 to better show the level of qualitative agreement. The quantitative discrepancies are shown in Fig. 8. Four points in the strain histories are marked out for the following discussion.

Because *De* is small, the transient response results mainly from kinematic hardening and thixotropy. Kinematic hardening leads to the stress transients within around three strain units after the shearing direction changes, as seen in Fig. 7(a). The thixotropic effects can be seen in Fig. 7(b): From the point a to b, the slope increases as the shear rate decreases to zero, which indicates a build-up of the viscosity; and the slope decreases again for the strain histories from the point b to c, as the shear rate reaches its negative maximum. A hysteresis loop is formed because the evolution of internal structures lags behind the change of the shear rate. The ML and MDT models lack the component of kinematic hardening. They, therefore, yield poor agreement with the experimental data in Fig. 7(b). The original IKH model includes kinematic hardening but it fails to quantitatively capture the thixotropic response due to the assumption of a constant plastic viscosity. In comparison, the IKH-V and ML-IKH models capture both kinematic hardening (through the evolution of *A*) and thixotropy (through the response of *λ*) and therefore give better predictions than the other three models.

The model predictions of the LAOS tests^{1} for additional test conditions are presented in Fig. 8. The IKH-V and ML-IKH models have very similar predictions in the LAOS tests although they show distinct results in step tests. MDT, IKH-V, and ML-IKH models all have accurate predictions when $\omega \gamma 0$ is large (in the ideal thixotropy region) but fail to predict experimental results, even qualitatively, when $\gamma 0=0.1$. The models fail because they oversimplify the complex viscoelastic response. The test fluid shows complex nonlinear viscoelastic response with multiple time scales while the models assume linear viscoelasticity with a single viscoelastic relaxation time. As $\gamma 0$ increases, the model predictions improve. In the tests for $\gamma 0=1$ and 10, the kinematic hardening effect is significant. The IKH-V and ML-IKH models capture this response and show better agreement with experimental results than the MDT model does.

### F. Summary

Table VI assesses the predictive ability of the MDT, ML, IKH-V, and ML-IKH models for transient shear rheology of thixotropic fluids. The ML-IKH model has four key aspects that result in the quantitative predictions of the complex TEVP rheology. They are (1) a spectrum of thixotropic structure parameters (corresponding model parameter: *β*), (2) nonlinear kinetic equations for thixotropic evolution (*a*, *n*), (3) the Armstrong-Frederick kinematic hardening rule ($kh$, *q, m*), and (4) the inclusion of viscoelasticity (*G*). Feature 1 generates stretched-exponential thixotropic evolution in step tests [Figs. 3(b), 4, and 11] and is critical in capturing the nonmonotonic stress transients in intermittent step tests (Fig. 5). Feature 2 enables the model to quantitatively predict step tests with various test conditions (Fig. 11). Feature 3 is critical in predicting the transient response induced by the change of shearing direction (Figs. 6–8). Feature 4 is important in capturing the fast transient response resulting from viscoelasticity [Figs. 3(b), 6, and 8]. In applications where only some of the features are of interest or importance, the model can be conveniently simplified by setting corresponding model parameters to specific values, as summarized in Table II.

. | MDT . | ML . | IKH-V . | ML-IKH . | Remarks . |
---|---|---|---|---|---|

Number of parameters | 13 | 9 | 9 | 12 | The original IKH model is modified by adding one model parameter, as explained in the beginning of Sec. V. |

Steady state flow curve | E | E | E | E | |

Step tests with small shear rate or stress | M | M | M | G | MDT fails to predict the initial elastic response. |

Step tests with large shear rate or stress | M | G | M | G | IKH-V and MDT fail to predict the stretched-exponential thixotropic relaxation. |

Intermittent step tests | P | M | P | G | IKH-V and MDT give poor predictions for intermediate $\Delta t$ |

Flow reversal tests | P | P | G | G | The IKH-V and ML-IKH models show good predictions for small shear rates. |

LAOS tests (large $\omega \gamma 0$) | M | M | G | G | The IKH-V and ML-IKH models have very similar predictions. |

LAOS tests (small $\omega \gamma 0$) | P | P | P | P | For the same $\omega \gamma 0$ value, model predictions improve as $\gamma 0$ increases. |

. | MDT . | ML . | IKH-V . | ML-IKH . | Remarks . |
---|---|---|---|---|---|

Number of parameters | 13 | 9 | 9 | 12 | The original IKH model is modified by adding one model parameter, as explained in the beginning of Sec. V. |

Steady state flow curve | E | E | E | E | |

Step tests with small shear rate or stress | M | M | M | G | MDT fails to predict the initial elastic response. |

Step tests with large shear rate or stress | M | G | M | G | IKH-V and MDT fail to predict the stretched-exponential thixotropic relaxation. |

Intermittent step tests | P | M | P | G | IKH-V and MDT give poor predictions for intermediate $\Delta t$ |

Flow reversal tests | P | P | G | G | The IKH-V and ML-IKH models show good predictions for small shear rates. |

LAOS tests (large $\omega \gamma 0$) | M | M | G | G | The IKH-V and ML-IKH models have very similar predictions. |

LAOS tests (small $\omega \gamma 0$) | P | P | P | P | For the same $\omega \gamma 0$ value, model predictions improve as $\gamma 0$ increases. |

## VI. PREDICTIONS OF THE TENSORIAL MODEL

### A. Comparison with the scalar model

Here, we demonstrate that the proposed tensorial model can reproduce the quantitative predictions of the scalar model. According to Eq. (31), the viscoelastic relaxation time is $\lambda \eta thi/G$. In simple shear flows, an effective Weissenberg number (*Wi*_{eff}) can be defined as

For infinitesimal *Wi*_{eff}, $Dp\u2248D$, and the flow parameter, *ϕ*, and consequently *λ*, are equal to those in the scalar model. Our model allows *Wi*_{eff} to remain small even for large shear rates because as shear rate increases, *λ* decreases. With some parameter values, e.g., *n* = 1 and *a* > 1.5 in the RC form, *Wi*_{eff} approaches zero for infinitely large $\gamma \u0307$.

With the same parameter values and finite *Wi*_{eff}, the shear stress and back strain predicted by the tensorial model are not identical to those of the scalar model. Consider the case where *θ* = 0, for which Eqs. (31) and (32) reduce to

In simple shear flows, the components of $\tau eff$ and $A$ are

$\sigma \xaf$ is related to $\gamma \u0307$ through

Therefore,

where $\kappa back,12=khA12(1+2A11)$. The evolution of $A11$ and $A12$ follows as:

In comparison, for the scalar model, $A\u0307=\gamma \u0307\u2212q\gamma \u0307A$ and

With different values of *k*_{h} and *q*, the tensorial evolution equation for **A** can reproduce the predictions of its scalar equivalent. To show this, Fig. 9 plots the predictions of the tensorial and scalar models in shear-start-up tests, with model parameters for both taken from Table IV for the SC form, except that for the tensorial model we take *k*_{h} = 0.3 Pa and *q* = 3. The *k*_{h} and *q* are determined so that $\kappa back,12$ and $khA$ have close results in flow histories with small *Wi*.

Figure 9 shows that as the shear rate increases, $\kappa back,12$ deviates from $khA$ and exhibits a maximum that increases with increasing shear rate. This discrepancy occurs because as shear rate and therefore *Wi*_{eff} increase, Eq. (43) is no longer valid and the first term in Eq. (31) is significant. Nevertheless, the stress difference resulting from this difference between $\kappa back,12$ and $khA$ remains insignificant relative to the viscous stress (for example, in test iv, the difference between $\kappa back,12$ and $khA$ is less than 0.3 Pa while the total shear stress is about 50 Pa). Consequently, the tensorial and scalar models give very similar predictions of shear stress and *λ*, as shown in Figs. 9(a) and 9(b). The tensorial model, therefore, can reproduce the scalar model's quantitative predictions in the test results shown in Figs. 3–8.

### B. Comparison with the tensorial IKH model

Here, we briefly compare our tensorial ML-IKH model with the tensorial IKH model, given in the supporting information of [15]. The tensorial IKH model is a 3D generalization of its scalar form based on the multiplicative decomposition of the deformation gradient tensor. It involves stress and strain measures defined in a “reference configuration” and an “intermediate configuration.” Such formulation belongs to the Lagrangian form and is more suitable for predicting the solid mechanics of TEVP materials. The tensorial ML-IKH model is based on the additive decomposition of the rate of deformation tensor. It does not depend on a reference configuration. The Cauchy stress tensor, its rate of change, the rate of deformation tensor, and the kinematic hardening variable can be written in a simple functional form as discussed in Sec. III B. This formulation belongs to the Eulerian form and is more suitable for most fluid flow simulations of TEVP fluids. The tensorial IKH model is in principle applicable for finite and large elastic deformations while the tensorial ML-IKH model assumes infinitesimal elastic deformation. This assumption is valid for weakly aggregated suspensions and it greatly simplifies the constitutive equations. With this assumption, the tensorial IKH model can also be written in Eulerian form. But it still has a drawback in that it predicts an unphysical oscillatory shear stress in simple shear flows for certain values of model parameters and shear rates. The oscillation results from its kinematic hardening rule, adopted from the theoretical work of Henann and Anand [52]. This disadvantage is however absent in the ML-IKH model.

For simplicity, consider the limit of zero yield stress and small elastic deformation. The tensorial IKH model in this limit can be written in a simple Eulerian form, as summarized in Table VII. In the limit of $\theta =0$, the two models further reduce to rigid, i.e., nonelastic, plastic models.

Additional limits . | ML-IKH . | IKH model . |
---|---|---|

None | $\tau =\tau s+2\eta mD\lambda \theta \u2207\u2009\u2009\tau s+\tau eff=2\lambda \eta thiD\tau eff=\tau \u2212\kappa back=2\eta thiDp\kappa back=kh(A+2A2)\omicron A=A\u22c5Dp+Dp\u22c5A+Dp\u2212qdpA$ | $\theta \omicron \tau +\tau eff=2\eta D\tau eff=\tau \u2212\tau back=2\eta Dp\tau back=khA=khlnA\u0303\omicron A\u0303=A\u0303\u22c5Dp+Dp\u22c5A\u0303\u2212qdpA\u0303\u22c5ln(A\u0303)$ |

Rigid plasticity $(\theta =0)$ | $\tau =2(\lambda \eta thi+\eta m)D+\kappa back\u2207A=D\u2212qdA$ | $\tau =2\eta D+\tau back\u2207A\u0303=\u2212qdA\u0303\u22c5ln(A\u0303)$ |

Remarks: | A is symmetric. A = 0 in the un-deformed reference configuration | $A\u0303$ is symmetric and unimodular. $A=lnA\u0303$ is symmetric and deviatoric. $A\u0303=1$ and $A=0$ in the undeformed reference configuration. |

Additional limits . | ML-IKH . | IKH model . |
---|---|---|

None | $\tau =\tau s+2\eta mD\lambda \theta \u2207\u2009\u2009\tau s+\tau eff=2\lambda \eta thiD\tau eff=\tau \u2212\kappa back=2\eta thiDp\kappa back=kh(A+2A2)\omicron A=A\u22c5Dp+Dp\u22c5A+Dp\u2212qdpA$ | $\theta \omicron \tau +\tau eff=2\eta D\tau eff=\tau \u2212\tau back=2\eta Dp\tau back=khA=khlnA\u0303\omicron A\u0303=A\u0303\u22c5Dp+Dp\u22c5A\u0303\u2212qdpA\u0303\u22c5ln(A\u0303)$ |

Rigid plasticity $(\theta =0)$ | $\tau =2(\lambda \eta thi+\eta m)D+\kappa back\u2207A=D\u2212qdA$ | $\tau =2\eta D+\tau back\u2207A\u0303=\u2212qdA\u0303\u22c5ln(A\u0303)$ |

Remarks: | A is symmetric. A = 0 in the un-deformed reference configuration | $A\u0303$ is symmetric and unimodular. $A=lnA\u0303$ is symmetric and deviatoric. $A\u0303=1$ and $A=0$ in the undeformed reference configuration. |

Figure 10 plots the predictions of $2A12$ by our model and its equivalent, $A12$, by the IKH model in the limit of rigid plasticity in shear-start-up tests. For small *q* values, the IKH model predicts a damped oscillation of $A12$ and therefore of the shear stress $(\eta \gamma \u0307+khA12)$. As *q* increases, this oscillation vanishes and the two models show the same predictions. The oscillation of shear stress is significant when both *q* and $\eta \gamma \u0307/kh$ are small. Note that *q* = 1 is used in this study to model the rheology of the fumed silica suspension, and it is yet unclear what the range of values that *q* takes in experimental TEVP systems.

## VII. CONCLUSION

The primary goal of this work is to develop a comprehensive constitutive model for the transient rheology of TEVP fluids. Such TEVP fluids show the following important rheological features: (1) Viscoelasticity and stretched-exponential thixotropic relaxation in step rate or step stress tests, (2) a nonmonotonic thixotropic response in intermittent flows, (3) kinematic hardening, and (4) linear viscoelasticity before yielding. We proposed here a phenomenological model that captures those rheological features by combining four key components that have not previously been unified in a single model: (1) Multiple thixotropic structure parameters, (2) nonlinear thixotropic kinetic equations, (3) the Armstrong-Frederick kinematic hardening rule, and (4) linear viscoelasticity. The proposed model is a generalization of our previous ML ideal thixotropic model and merger of this with the IKH model. This new model has 12 parameters, four more than the IKH model, where the additional parameters introduce multiple thixotropic structural parameters which allow the structure to relax in a stretched exponential form, as well nonlinearities into the evolution of structure, and leading to more accurate predictions in step rate or step stress tests. Such a large number of parameters is necessary to accurately capture the TEVP rheology and is common in thixotropic models [15,24,26,27,31,58].

For limiting values of model parameters, the ML-IKH model reduces to several rheological models, such as the Oldroyd-B model (when thixotropy and plasticity are absent), the Bingham model (when viscoelasticity, thixotropy, and back stress are absent), Saramito's elastoviscoplastic model [51] (when thixotropy and back stress are absent), Moore's ideal thixotropic model [71] (when viscoelasticity and back stress are absent), and more examples are summarized in Table II. This flexibility allows the users to conveniently simplify it as needed for desired applications.

We comparatively examined this model over a wide range of test conditions for steady shear, step shear rate, step stress, intermittent shear, shear reversal, and LAOS experiments. The experimental data are from a well-studied thixotropic fumed silica suspension. The results of the comparative testing are summarized in Table VI. This comparative study demonstrates the success of the four key model components in capturing several typical rheological features of TEVP materials. Multiple thixotropic structure parameters and their nonlinear kinetic equations are critical for predicting the rheology in step tests and intermittent shear tests. The Armstrong-Frederick kinematic hardening rule is important in predicting the flow reversal tests. The model components for linear viscoelasticity capture the fast elastic response in step tests.

We also presented a tensorial form of this model based on the additive decomposition of the rate of deformation tensor. It is written in the Eulerian form and suitable for the fluid flow simulation of thixotropic fluids. We demonstrated that this tensorial model can reproduce the quantitative predictions of the scalar model and does not show unphysical oscillation of shear stress in simple shearing flows. It is frame-invariant and satisfies the second law of thermodynamics (see Appendix B for details).

This model has following limitations. First, the ML approach accounting for thixotropy assumes a constant stretching factor (*β*) while the experiments in [37] have shown that *β* varies as test condition changes. Breakdown experiments with large final shear rates exhibit small *β* values (e.g., 0.3); conversely, build-up tests with small final shear rates give large *β* values (e.g., 1.5). This limitation may lead to poor predictions of the thixotropic evolution when the shear rate is small. Second, this model fail to predict the delayed yielding behavior that some TEVP materials exhibit. Third, this model assumes linear viscoelasticity with a single time scale for viscoelastic relaxation. The real TEVP fluids, however, often exhibit complex nonlinear viscoelasticity. Fourth, the model assumes that the yield stress depends on structural anisotropy (through the Armstrong-Frederick kinematic hardening rule of back stress) but the viscosity does not, which needs to be further examined. Furthermore, the tensorial form of this model in principle can be used in arbitrary 3D flow conditions. But it has hitherto been tested only in simple shear flows because of the lack of nonshearing rheological experiments on TEVP materials. We leave it for further work to overcome these limitations.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge the support from the Procter & Gamble Company. M.J.S. acknowledges partial support from NSF CBET 1232937 and R.G.L. from NSF CBET 0853662. 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 National Science Foundation (NSF).

## NOMENCLATURE

*a*order of flow parameter in break-down term

*A*kinematic hardening structure parameter

**A**3D generalization of

*A**b*order of flow parameter in shear-induced build-up term

**B**_{e}elastic Finger tensor

*C*_{i}coefficient of

*λ*_{i}**D**rate of deformation tensor

*D*_{i}prefactor of

*λ*in its kinetic equation_{i}**D**_{e}**, D**_{p}elastic and plastic contribution of

**D****F**deformation gradient

**F**_{e},**F**_{p}elastic and plastic components of

**F***G*elastic modulus

*k*_{1}kinetic constant of break-down term

*k*_{2}kinetic constant of shear-induced build-up term

*k*_{3}kinetic constant of Brownian build-up term

*k*_{h}back stress coefficient

*k*_{y}yield stress coefficient

*K*_{1}*k*_{1}/*k*_{3}*K*_{2}*k*_{2}/*k*_{3}- $\gamma \u0307$
shear rate

*m*order of

*λ*in break-down term*n*order of

*λ*in build-up terms*N*number of

*λ*_{i}*t*time

**W**vorticity tensor

*β*stretching factor

- $\u2207v$
velocity gradient

*η*_{m}limiting high shear rate viscosity

*η*_{thi}thixotropic viscosity increment

*ϕ*flow parameter

- $\gamma \u0307$
shear rate

*λ*_{i}substructure parameter

*λ*structure parameter

- ψ
free energy

- ψ
_{e}, ψ_{p} elastic and plastic free energies

*σ*stress

*σ*_{∞}equilibrium value of stress

*σ*_{o}initial stress

*σ*_{t}final stress

*σ*_{y0}apparent yield stress

**τ**extra stress tensor

**τ**_{m},**τ**_{s}stress contribution from medium and internal structures

**τ**_{eff}effective stress tensor

**τ**_{back}back stress tensor

**κ**_{back}intermediate back stress tensor

*Δt*duration of the intermediate step

### APPENDIX A: MODEL PREDICTIONS IN ADDITIONAL STEP TESTS

Here, we test the models in step tests over a range of test conditions. Figures 11(a) and 11(b) plot the model predictions in nine step-shear-rate tests. The RC and SC forms of the ML-IKH model give similar predictions that are both in quantitative agreement with the experimental data, as shown in Fig. 11(a). The MDT model predicts those tests well but tends to predict narrower stress transients than the experimental results. The IKH model provides only qualitatively correct predictions.

### APPENDIX B: FREE ENERGY IMBALANCE

In this section, we demonstrate that the tensorial ML-IKH model satisfies the free-energy imbalance, i.e., obeys the second law of thermodynamics. The total free energy ($\Psi $) has the following form:

where the elastic energy is described the incompressible neo-Hookean model. Here the first and second terms on the right side are, respectively, the elastic ($\Psi e$) and plastic ($\Psi p$) free energies. $Be$ denotes the elastic finger tensor, defined as $Be=Fe\u22c5FeT$, where $Fe$ is the elastic component of the deformation gradient $F$ according to the Kroner decomposition [67], $F=Fe\u22c5Fp$. Here, $F$, $Fe$, and $Fp$ are defined as

where **X**, $x\xaf$, and **x** are, respectively, position vectors in the reference, intermediate, and current configurations. The velocity gradient tensor ($L=(\u2207v)T$) as well as its elastic and plastic components are, respectively, $L=F\u0307\u22c5F\u22121$, $Le=F\u0307e\u22c5Fe\u22121$, and $Lp=Fe\u22c5F\u0307p\u22c5Fp\u22121\u22c5Fe\u22121$.

For incompressible flows, the free-energy imbalance requires that $\tau :D\u2212\Psi \u0307>0$ is always valid, where $\tau =\tau m+\tau s$. The rate of change of free energy is $\Psi \u0307=\Psi \u0307e+\Psi \u0307p$ where $\Psi \u0307e=G2\delta :B\u0307e$, and $\Psi \u0307p=khA:A\u0307$. We expand $B\u0307e$ as follows:

Therefore,

which gives

We next derive the expression for $\u2009\kappa back:Dp$, which will be used in Eq. (B7)

According to Eq. (29), $\tau eff:Dp$ is nonnegative. And $\tau back:A=khA:A\u22650$. Therefore, the free energy imbalance is always satisfied.

^{1}

The ML model works well only when $\omega \gamma 0$ is large ($\u2248100\u2009s\u22121$ in this study) so that the viscoelastic and plastic response are insignificant. For brevity, we do not include the ML model in Fig. 8. Its predictions are provided in Fig. S6.