3.1 Fabric tensor

In a non-Brownian hard-sphere suspension, particles interact solely through hydrodynamic and hard contact forces. While contact forces are short-range pairwise interactions, hydrodynamic forces result from the superposition of algebraically decaying solvent velocity fields created by the particles’ motion, and as such have a long-range and many-body nature. In a concentrated system, however, most of the hydrodynamic tractions on the particle surfaces come from lubrication flows within the narrow interparticle gaps (Frankel & Acrivos Reference Frankel and Acrivos1967). These lubrication forces are pairwise and short-ranged in nature. Therefore, in the dense regime, most of the stress comes from pairwise short-range interactions (Ball & Melrose Reference Ball and Melrose1997; Mari et al. Reference Mari, Seto, Morris and Denn2015); indeed, the long-range hydrodynamics is omitted from our simulations, as described above. The pairwise nature of dominant forces will give a special importance for stress prediction to the pair correlation function $g(\boldsymbol{r})$ , defined as the conditional probability to find a particle whose centre is at $\boldsymbol{r}$ knowing that there is a particle at $\mathbf{0}$ , normalized by the average number density $\unicode[STIX]{x1D70C}$ . Moreover, because the dominant forces are also short-ranged and $g(\boldsymbol{r})$ is peaked at or very near contact (Morris & Katyal Reference Morris and Katyal2002; Nazockdast & Morris Reference Nazockdast and Morris2012, Reference Nazockdast and Morris2013; Mari et al. Reference Mari, Seto, Morris and Denn2014), we can expect that only the near-contact part of $g(\boldsymbol{r})$ will play a substantial role. Hence we can consider the distribution of near-contact orientations $P(\boldsymbol{p})$ (with $\boldsymbol{p}\in S^{2}$ ), which keeps only the orientational information of interactions that determine the stress state of the suspension.

Building a time-evolution dynamics for $P(\boldsymbol{p})$ (let alone $g(\boldsymbol{r})$ ) is, however, very challenging for a non-Brownian suspension, and quite possibly over-reaching in the context of building a constitutive model. Indeed, the stress $\unicode[STIX]{x1D72E}$ is a symmetric second-rank tensor with six independent components, and in consequence contains far less information than the entire $P(\boldsymbol{p})$ . Therefore, rather than attempt to model the full $P(\boldsymbol{p})$ (or $g(\boldsymbol{r})$ (Nazockdast & Morris Reference Nazockdast and Morris2012, Reference Nazockdast and Morris2013)) directly, previous attempts model the symmetric and traceless fabric tensor $\unicode[STIX]{x1D64C}:=\langle \boldsymbol{p}\boldsymbol{p}\rangle -{\textstyle \frac{1}{3}}\unicode[STIX]{x1D644}$ . This is the deviatoric second moment of $P(\boldsymbol{p})$ , and the lowest-order moment carrying non-trivial structural information. It also corresponds to the second-order term in the Laplace spherical harmonic expansion of $P(\boldsymbol{p})$ in tensor form (Kanatani Reference Kanatani1984):

(3.1) $$\begin{eqnarray}\displaystyle P(\boldsymbol{p})\approx \frac{1}{4\unicode[STIX]{x03C0}}\left(1+\frac{15}{2}\unicode[STIX]{x1D64C}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}\right). & & \displaystyle\end{eqnarray}$$

[For a symmetric distribution on the unit sphere like $P(\boldsymbol{p})$ , it is possible to define a tensorial expansion in powers of $\boldsymbol{p}$ as $P(\boldsymbol{p})=(1/(4\unicode[STIX]{x03C0}))\sum _{s=0}^{\infty }\unicode[STIX]{x1D60B}_{i_{1},\ldots ,i_{2s}}p_{i_{1}}\cdots p_{i_{2s}}$ such that the term of order $n=2s$ in the tensorial expansion is the projection of $P(\boldsymbol{p})$ on the $(2n+1)$ -dimensional space of Laplace spherical harmonics of order $n$ (Kanatani Reference Kanatani1984). Put differently, the order- $n$ term of the tensorial expansion is the sum of the terms of order $n$ in the Laplace spherical harmonics expansion. The tensors $\unicode[STIX]{x1D60B}_{i_{1},\ldots ,i_{n}}$ are deviatoric, hence contain $2n+1$ independent components, and are mutually independent, a consequence of the orthogonality of the Laplace spherical harmonics basis. The first-term tensor (zeroth rank) of the tensorial expansion is 1 (which enforces the normalization of $P(\boldsymbol{p})$ ), the tensor of the second term is ${\textstyle \frac{15}{2}}\unicode[STIX]{x1D64C}$ .]

Retaining the fabric tensor $\unicode[STIX]{x1D64C}$ as a proxy for the full microstructure is indeed the standard choice made in previous attempts to build a constitutive model of suspensions (see Hinch & Leal Reference Hinch and Leal1975; Phan-Thien Reference Phan-Thien1995; Phan-Thien, Fan & Khoo Reference Phan-Thien, Fan and Khoo1999; Goddard Reference Goddard2006; Stickel et al. Reference Stickel, Phillips and Powell2006). Note that there is also a long history of representing the microstructure via a second-rank tensor in polymeric systems (Larson Reference Larson2013) and in dry granular materials (Magnanimo & Luding Reference Magnanimo and Luding2011; Sun & Sundaresan Reference Sun and Sundaresan2011; Goddard Reference Goddard2014).

3.2 Dynamics: Hand equation

Once one decides to write the dynamics of the microstructure as a closed set of ordinary differential equations (ODEs) for the time evolution of the fabric tensor, symmetries and frame indifference (also known as material objectivity) constrain quite strongly the functional forms involved. In the absence of inertia, systems are invariant under time-dependent translations and rigid rotations, a property known as frame indifference. Using a representation theorem and enforcing frame indifference, one can write a general evolution equation for a symmetric $3\times 3$ tensor (Hand Reference Hand1962). Specifically, if $\text{d}\unicode[STIX]{x1D64C}/\text{d}t$ depends only on $\unicode[STIX]{x1D64C}$ itself and on the rate-of-strain tensor $\unicode[STIX]{x1D646}:=(\unicode[STIX]{x2202}_{j}v_{i})$ , where $\boldsymbol{v}$ is the local suspension velocity, we have

(3.2) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D64C}}{\text{d}t} & = & \displaystyle \unicode[STIX]{x1D652}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}-\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\unicode[STIX]{x1D652}+\tilde{\unicode[STIX]{x1D6FC}}_{0}\unicode[STIX]{x1D644}+\tilde{\unicode[STIX]{x1D6FC}}_{1}\unicode[STIX]{x1D64C}+\tilde{\unicode[STIX]{x1D6FC}}_{2}\unicode[STIX]{x1D640}+\tilde{\unicode[STIX]{x1D6FC}}_{3}\unicode[STIX]{x1D64C}^{2}+\tilde{\unicode[STIX]{x1D6FC}}_{4}\unicode[STIX]{x1D640}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\tilde{\unicode[STIX]{x1D6FC}}_{5}(\unicode[STIX]{x1D640}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\unicode[STIX]{x1D640})+\tilde{\unicode[STIX]{x1D6FC}}_{6}(\unicode[STIX]{x1D640}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}^{2}+\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D640})\nonumber\\ \displaystyle & & \displaystyle +\,\tilde{\unicode[STIX]{x1D6FC}}_{7}(\unicode[STIX]{x1D640}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{2})+\tilde{\unicode[STIX]{x1D6FC}}_{8}(\unicode[STIX]{x1D640}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}^{2}+\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{2}),\end{eqnarray}$$

where $\unicode[STIX]{x1D652}$ and $\unicode[STIX]{x1D640}$ are respectively the antisymmetric and symmetric parts of $\unicode[STIX]{x1D646}$ , and where the scalar coefficients $\tilde{\unicode[STIX]{x1D6FC}}_{i}$ are analytic functions of the invariants $I_{1}:=\text{Tr}[\unicode[STIX]{x1D64C}]$ , $I_{2}:=\text{Tr}[\unicode[STIX]{x1D64C}^{2}]$ , $I_{3}:=\text{Tr}[\unicode[STIX]{x1D64C}^{3}]$ , $I_{4}:=\text{Tr}[\unicode[STIX]{x1D640}]$ , $I_{5}:=\text{Tr}[\unicode[STIX]{x1D640}^{2}]$ , $I_{6}:=\text{Tr}[\unicode[STIX]{x1D640}^{3}]$ , $I_{7}:=\text{Tr}[\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\unicode[STIX]{x1D640}]$ , $I_{8}:=\text{Tr}[\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D640}]$ , $I_{9}:=\text{Tr}[\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{2}]$ and $I_{10}:=\text{Tr}[\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D640}^{2}]$ . Note that (3.2) follows only from frame indifference and the fact that $\unicode[STIX]{x1D64C}$ is a symmetric $3\times 3$ tensor. It does not rely on any specific physical property of the system nor on an expansion to some order in $\unicode[STIX]{x1D64C}$ or $\unicode[STIX]{x1D640}$ . As such, any closed ODE for $\text{d}\unicode[STIX]{x1D64C}/\text{d}t$ analytic in $\unicode[STIX]{x1D64C}$ and $\unicode[STIX]{x1D646}$ has to be of the form of (3.2).

In our case, $\unicode[STIX]{x1D64C}$ is traceless, which implies that $I_{1}=0$ and $\tilde{\unicode[STIX]{x1D6FC}}_{0}$ is prescribed by the value of the other coefficients $\tilde{\unicode[STIX]{x1D6FC}}_{i}$ . Incompressibility sets $\text{Tr}[\unicode[STIX]{x1D640}]=:I_{4}$ to zero too. We must, however, take special care when adopting Hand’s result for a non-Brownian hard-particle suspension. For a system with rate-independent particle trajectories, the fabric evolution must be proportional to the absolute value of the strain rate $|\dot{\unicode[STIX]{x1D6FE}}|=\sqrt{2\unicode[STIX]{x1D640}\boldsymbol{ :}\unicode[STIX]{x1D640}}$ . This is both necessary, because $\dot{\unicode[STIX]{x1D6FE}}$ is the only time scale, and permissible, because $\dot{\unicode[STIX]{x1D6FE}}=0$ corresponds to a singular case in which nothing evolves at all, so that analyticity must be defined on a domain in which $\dot{\unicode[STIX]{x1D6FE}}\neq 0$ . This is why the basic character of rate-independent rheology involves a linear dependence on $\dot{\unicode[STIX]{x1D6FE}}$ (for shear stress) or $|\dot{\unicode[STIX]{x1D6FE}}|$ (for normal stresses). The fabric evolution is then analytic in $\unicode[STIX]{x1D64C}$ and $\unicode[STIX]{x1D646}$ , except at $\unicode[STIX]{x1D646}=0$ (or equivalently $\dot{\unicode[STIX]{x1D6FE}}=0$ ). Also, at least for the simple shear flows studied below, we can change variable from time to the accumulated strain $\unicode[STIX]{x1D6FE}$ . Hence, by introducing $\hat{\unicode[STIX]{x1D640}}=\unicode[STIX]{x1D640}/|\dot{\unicode[STIX]{x1D6FE}}|$ , $\hat{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}/|\dot{\unicode[STIX]{x1D6FE}}|$ and $\dot{\unicode[STIX]{x1D64C}}=\text{d}\unicode[STIX]{x1D64C}/\text{d}\unicode[STIX]{x1D6FE}$ , we introduce a rate-independent Hand-like equation as

(3.3) $$\begin{eqnarray}\displaystyle \operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})\dot{\unicode[STIX]{x1D64C}} & = & \displaystyle \hat{\unicode[STIX]{x1D652}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}-\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D652}}+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D6FC}_{2}\hat{\unicode[STIX]{x1D640}}+\unicode[STIX]{x1D6FC}_{3}\unicode[STIX]{x1D64C}^{2}+\unicode[STIX]{x1D6FC}_{4}\hat{\unicode[STIX]{x1D640}}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{5}(\hat{\unicode[STIX]{x1D640}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}})+\unicode[STIX]{x1D6FC}_{6}(\hat{\unicode[STIX]{x1D640}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}^{2}+\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}})\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{7}(\hat{\unicode[STIX]{x1D640}}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}}^{2})+\unicode[STIX]{x1D6FC}_{8}(\hat{\unicode[STIX]{x1D640}}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}^{2}+\unicode[STIX]{x1D64C}^{2}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}}^{2})\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D6FC}_{3}\hat{I} _{2}+\unicode[STIX]{x1D6FC}_{4}\hat{I} _{5}+2\unicode[STIX]{x1D6FC}_{5}\hat{I} _{7}+2\unicode[STIX]{x1D6FC}_{6}\hat{I} _{8}+2\unicode[STIX]{x1D6FC}_{7}\hat{I} _{9}+2\unicode[STIX]{x1D6FC}_{8}\hat{I} _{10})\unicode[STIX]{x1D644},\end{eqnarray}$$

where the $\unicode[STIX]{x1D6FC}_{i}$ are now dimensionless, analytic functions of the nine invariants $\hat{I} _{k}$ involving $\unicode[STIX]{x1D64C}$ and $\hat{\unicode[STIX]{x1D640}}$ .

Hand’s equation in the form (3.3) provides the boundaries within which one can build a dynamics for $\unicode[STIX]{x1D64C}$ . It is nonetheless not practically useful as it stands, as there is no prescription for the coefficients $\unicode[STIX]{x1D6FC}_{i}$ , and the space of allowed dynamical equations is still infinite-dimensional. At this point, one usually has to introduce further assumptions about the specific nature of the dynamics in order to reduce the number of free parameters to a handful or so, based on physical motivations or simply tractability of the model.

4.1 A family of fabric tensors

There is some freedom in the exact definition of $\unicode[STIX]{x1D64C}$ . We would like $\unicode[STIX]{x1D64C}$ to be an average over the interactions that are relevant to the suspension stress, which we know are short-range. What do we mean by short-range interaction precisely? One could consider a $\unicode[STIX]{x1D64C}$ including only strict contact interactions, or including every pair of particles sharing an edge in a Delaunay tessellation (i.e. particles occupying neighbouring cells in a Voronoi tessellation), or including every pair of particles separated by a gap smaller than $\unicode[STIX]{x1D716}$ . We adopt the latter choice and define the set of directions of interactions closer than a gap $\unicode[STIX]{x1D716}$ as $\unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D716}}=\{\,\boldsymbol{p}_{ij}~\forall \,i,j\mid 2(r_{ij}/(a_{i}+a_{j})-1)<\unicode[STIX]{x1D716}\}$ , where $i$ and $j$ are particle labels, $r\boldsymbol{p}$ is a centre-to-centre vector and $a$ is a particle radius. We then obtain a family of traceless fabric tensors $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ as

(4.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}=\langle \boldsymbol{p}\boldsymbol{p}\rangle _{\boldsymbol{ p}\in \unicode[STIX]{x1D6E4}^{\unicode[STIX]{x1D716}}}-{\textstyle \frac{1}{3}}\unicode[STIX]{x1D644}. & & \displaystyle\end{eqnarray}$$

Using numerical simulations, we can study the dynamics of $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ for any $\unicode[STIX]{x1D716}$ . In figure 1, we show that the dynamics of $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}=0}$ (i.e. including only strict contacts, which correspond to interactions with a small negative gap thanks to the finite deformability of the particles) and $\unicode[STIX]{x1D72E}$ share several features upon reversal. In particular, they both show a discontinuity exactly at reversal, followed by a relaxation back to steady state over the same strain scale, here roughly 2 strain units for $\unicode[STIX]{x1D719}=0.5$ . (Strictly speaking, because of the slight deformability of the particles, the apparent discontinuity in $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}=0}$ at reversal is actually a fast relaxation taking place during the first 1 % of strain after reversal; see figure 1(b) inset.) Because of the symmetry of the simple shear flow, both admit non-zero values only for their diagonal components and for their shear components in the velocity-gradient plane, $Q_{12}$ and $\unicode[STIX]{x1D6F4}_{12}$ . The steady-state values of these components have the same parity under strain reversal, with sign reversal for shear and none for the diagonal elements, providing further motivation for the use of $\unicode[STIX]{x1D64C}$ as a proxy for the microstructure in a constitutive model for the stress $\unicode[STIX]{x1D72E}$ .

Figure 1.

Stress and fabric data from numerical simulations. (a) Stress tensor components $\unicode[STIX]{x1D6F4}_{12}$ , $N_{1}=\unicode[STIX]{x1D6F4}_{11}-\unicode[STIX]{x1D6F4}_{22}$ and $N_{2}=\unicode[STIX]{x1D6F4}_{22}-\unicode[STIX]{x1D6F4}_{33}$ divided by the absolute value of the shear rate $|\dot{\unicode[STIX]{x1D6FE}}|$ upon shear reversal at $\unicode[STIX]{x1D6FE}=0$ for $\unicode[STIX]{x1D719}=0.5$ . Components $\unicode[STIX]{x1D6F4}_{13}$ and $\unicode[STIX]{x1D6F4}_{23}$ are not shown, as they vanish in simple shear by symmetry. Thick dark-shaded lines are the averaged data, while the light shaded area around each curve is the standard deviation obtained from the individual shear reversals. All components show a discontinuity upon reversal, arising from lubrication and/or contact forces. (b) Fabric tensor components (with $\unicode[STIX]{x1D716}=0$ , describing full contacts only) for the same conditions, with averages in thick lines and standard deviation in shaded areas. Not shown are $Q_{13}=Q_{23}=0$ and $Q_{33}=-Q_{11}-Q_{22}$ . Insets: zoom-in on strains near the point of reversal.

Figure 2.

Fabric tensor components upon reversal at $\unicode[STIX]{x1D6FE}=0$ for $\unicode[STIX]{x1D719}=0.5$ and for three different cut-off lengths, (a)  $\unicode[STIX]{x1D716}=0$ , (b)  $\unicode[STIX]{x1D716}=0.02$ and (c)  $\unicode[STIX]{x1D716}=0.1$ , from numerical simulations. The fabric tensor built on contacts (i.e. with $\unicode[STIX]{x1D716}=0$  (a)) is discontinuous upon reversal, as a finite fraction of contacts disappear when the flow is reversed. However, when the cut-off length is finite (i.e. (b) and (c)), the discontinuity disappears, and instead there is a delay followed by a rapid evolution of the $Q_{ij}$ from their pre-reversal steady-state values to their maximum ( $Q_{12}$ and $Q_{22}$ ) or minimum ( $Q_{11}$ ) values, both phases occurring over a strain of order  $\unicode[STIX]{x1D716}$ .

However, despite the appeal of the similarity with the dynamics of the stress, $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}=0}$ is discontinuous at reversal whereas solutions of (3.3) are continuous. In figure 2, we show instead that, for $\unicode[STIX]{x1D716}>0$ , $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ is continuous upon reversal and thus possibly described by an evolution like (3.3), and otherwise keeps most of the qualitative features of $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}=0}$ . The initial singular reversal behaviour is smoothed over a finite strain scale which is increasing with $\unicode[STIX]{x1D716}$ and saturates for large $\unicode[STIX]{x1D716}$ to a strain $O(1)$ . Separately, the overall amplitude of $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ decreases with  $\unicode[STIX]{x1D716}$ .

An appealing strategy is then to assume that knowledge of the fabric of near-contacts, $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}>0}$ , along with knowledge of the flow itself, is enough to predict the evolution of the stress tensor $\unicode[STIX]{x1D72E}$ . That is, a continuous time evolution of $\unicode[STIX]{x1D64C}$ obeying an ODE in the Hand class might be married with a rule for constructing $\unicode[STIX]{x1D72E}(\unicode[STIX]{x1D64C},\unicode[STIX]{x1D646})$ that, through its dependence on $\unicode[STIX]{x1D646}$ , is explicitly alert to discontinuities at the moment of strain reversal. (This rule could build in information about the sudden opening of contacts on scales below $\unicode[STIX]{x1D716}$ , for instance.) An alternative strategy is to build a time evolution for the true contact fabric ( $\unicode[STIX]{x1D716}=0$ ) based on two relaxation times, one large corresponding to the time scale of reorganization of the microstructure, and one extremely small to describe the instantaneous changes in the contacts upon flow changes (Goddard Reference Goddard2006). However, with this approach, the fit to actual data imposes the parameters associated with the short time relaxation to take values whose physical interpretation is complex, including, for instance, transiently negative diffusion coefficients (Goddard Reference Goddard2008).

We stress that the cut-off length $\unicode[STIX]{x1D716}$ is not related to the lubrication regularization length $\unicode[STIX]{x1D6FF}$ nor to the maximum allowed overlap in the contact model. Moreover, it appears that, because $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ is getting smoothed out with increasing $\unicode[STIX]{x1D716}$ , the effects of contact overlaps and $\unicode[STIX]{x1D6FF}$ , which are only visible in particle trajectories at very early strains after reversal, will be gradually hidden in the time evolution of $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ for increasing $\unicode[STIX]{x1D716}$ . In practice, we find that $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ becomes insensitive to a doubling of $\unicode[STIX]{x1D6FF}$ for $\unicode[STIX]{x1D716}$ as small as 0.01. With this as the modelling strategy, picking the best cut-off length $\unicode[STIX]{x1D716}$ to build a constitutive model is then a compromise between the accuracy needed to admit a quantitative stress rule from $\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}}$ , and the smoothness required for an ODE closure in the $\unicode[STIX]{x1D64C}$ dynamics. In this work, we do not address the issue of the stress rule further, addressing solely the dynamics of $\unicode[STIX]{x1D64C}$ itself. We will focus our numerical tests of fabric evolution models on the case $\unicode[STIX]{x1D716}=0.02$ , adopting the default notation $\unicode[STIX]{x1D64C}:=\unicode[STIX]{x1D64C}^{\unicode[STIX]{x1D716}=0.02}$ . But, as shown in figures 2 and 3, the qualitative features of the dynamics are shared across all values of $\unicode[STIX]{x1D716}$ , so that most of what we will learn by trying to fit the case $\unicode[STIX]{x1D716}=0.02$ readily extends to other values of  $\unicode[STIX]{x1D716}$ .

4.2 Weak anisotropy: small- $\unicode[STIX]{x1D64C}$ models

Figure 3.

Fabric tensor components upon reversal at $\unicode[STIX]{x1D6FE}=0$ for three volume fractions, $\unicode[STIX]{x1D719}=0.4$ (a–c), 0.5 (d–f) and 0.55 (g–i), and for three different cut-off scales, $\unicode[STIX]{x1D716}=0.01$ (a,d,g), 0.02 (b,e,h) and 0.1 (c,f,i), from numerical simulations.

At large volume fractions, the steric constraints impose that every particle in the suspension is surrounded by close neighbours in any direction, with possibly only one direction (the compressional axis in simple shear) being slightly more crowded than the others. We can then expect that the anisotropy of the microstructure is decreasing with increasing  $\unicode[STIX]{x1D719}$ .

Figure 4.

Anisotropy $A$ as a function of post-reversal strain for three volume fractions, $\unicode[STIX]{x1D719}=0.4$ , 0.5 and 0.55, from numerical simulations.

As evidence for this, given a distribution of near-contacts $P(\boldsymbol{p})$ , we quantify its anisotropy by defining $A:=(4\unicode[STIX]{x03C0})^{-1}\int _{S^{2}}[4\unicode[STIX]{x03C0}P(\boldsymbol{p})-1]^{2}\,\text{d}\unicode[STIX]{x1D6FA}$ . Figure 4 shows the post-reversal evolution of $A$ for various volume fractions $\unicode[STIX]{x1D719}$ , and makes it clear that $A$ does indeed decrease with increasing $\unicode[STIX]{x1D719}$ . Similarly, in figure 3, we show the dependence of $\unicode[STIX]{x1D64C}$ on $\unicode[STIX]{x1D719}$ , for several values of $\unicode[STIX]{x1D716}$ , which demonstrates for $\unicode[STIX]{x1D64C}$ the corresponding lowering of the anisotropy as $\unicode[STIX]{x1D719}$ increases. To quantify further how small $\unicode[STIX]{x1D64C}$ actually is, we can look at a norm such as $|\unicode[STIX]{x1D64C}|=\sqrt{\text{Tr}[\unicode[STIX]{x1D64C}^{2}]}$ . It is easy to show that $|\unicode[STIX]{x1D64C}|\leqslant \sqrt{2/3}\approx 0.816$ , with equality corresponding to having all the $\boldsymbol{p}$ in the same direction, $P(\boldsymbol{p})=\unicode[STIX]{x1D6FF}(\boldsymbol{p}-\boldsymbol{p}_{0})$ . In contrast with this saturation value, we observe in steady state for $\unicode[STIX]{x1D719}=0.4$ and $\unicode[STIX]{x1D716}=0.01$ (figure 3 a), $|\unicode[STIX]{x1D64C}|\approx 0.048$ , and for $\unicode[STIX]{x1D719}=0.55$ and $\unicode[STIX]{x1D716}=0.01$ (figure 3 g),  $|\unicode[STIX]{x1D64C}|\approx 0.0022$ . These numbers, which confirm a weak anisotropy limit, become even smaller if we increase $\unicode[STIX]{x1D716}$ . A reasonable assumption, based on the smallness of these values, is to restrict the evolution equation (3.3) to low orders in $\unicode[STIX]{x1D64C}$ , allowing a Taylor expansion of the dynamics around the isotropic state. In what follows, we will examine the resulting dynamical equations at successively higher orders in  $\unicode[STIX]{x1D64C}$ .

4.3 Hand equation in simple shear

From now on, we restrict our discussion to a simple shear flow with $\hat{K}_{ij}=0$ except $\hat{K}_{12}=\pm 1/2$ , which is the simulation case. Because of the sparsity of $\hat{\unicode[STIX]{x1D646}}$ in this geometry, not all the components on the right-hand side of (3.3) are linearly independent, and some of the $\unicode[STIX]{x1D6FC}_{i}$ are redundant. Hence we recast the tensorial Hand equation for rate-independent systems (3.3) in a non-redundant component form (see appendix B) as

(4.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\dot{Q}_{-}=P_{-}[Q_{+},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},Q_{-}^{2}]\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{-}+2Q_{12},\\ \dot{Q}_{+}=P_{+}[Q_{+},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},Q_{-}^{2}]\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}}),\\ \dot{Q}_{12}=P_{12}[Q_{+},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},Q_{-}^{2}]-\frac{1}{2}Q_{-},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with functions $P_{-}$ , $P_{+}$ and $P_{12}$ analytic in their arguments, where $Q_{\pm }:=Q_{11}\pm Q_{22}$ .

We see from (4.2) that, when fitting these equations to simulation data, the basis of choice is $\{Q_{+},Q_{-},Q_{12}\}$ , rather than the naive basis $\{Q_{11},Q_{22},Q_{12}\}$ . Choosing the $\{Q_{+},Q_{-},Q_{12}\}$ basis when plotting fits has the effect of highlighting the important features of the data that a given model is trying to fit. For reversal protocols in particular, apparent qualitative agreement between model and data for $Q_{11}(\unicode[STIX]{x1D6FE})$ and $Q_{22}(\unicode[STIX]{x1D6FE})$ can be exposed as clear disagreement in this basis.

5.1 Linear Hand equation

We begin by assuming the fabric evolves via a frame-indifferent ODE linear in $\unicode[STIX]{x1D64C}$ . From (3.2), the most general linear model is

(5.1) $$\begin{eqnarray}\displaystyle \operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})\dot{\unicode[STIX]{x1D64C}} & = & \displaystyle \hat{\unicode[STIX]{x1D652}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}-\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D652}}+\unicode[STIX]{x1D6FC}_{1}\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D6FC}_{2}\hat{\unicode[STIX]{x1D640}}+\unicode[STIX]{x1D6FC}_{4}\hat{\unicode[STIX]{x1D640}}^{2}\nonumber\\ \displaystyle & & \displaystyle +\,\unicode[STIX]{x1D6FC}_{5}(\hat{\unicode[STIX]{x1D640}}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}})+\unicode[STIX]{x1D6FC}_{7}(\hat{\unicode[STIX]{x1D640}}^{2}\boldsymbol{\cdot }\unicode[STIX]{x1D64C}+\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}}^{2})\nonumber\\ \displaystyle & & \displaystyle -\,{\textstyle \frac{1}{3}}(\unicode[STIX]{x1D6FC}_{4}\hat{I} _{5}+2\unicode[STIX]{x1D6FC}_{5}\hat{I} _{7}+2\unicode[STIX]{x1D6FC}_{7}\hat{I} _{9})\unicode[STIX]{x1D644},\end{eqnarray}$$

with

(5.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FC}_{1}=x_{10},\\ \unicode[STIX]{x1D6FC}_{2}=x_{20}+x_{27}\hat{I} _{7}+x_{29}\hat{I} _{9},\\ \unicode[STIX]{x1D6FC}_{4}=x_{40}+x_{47}\hat{I} _{7}+x_{49}\hat{I} _{9},\\ \unicode[STIX]{x1D6FC}_{5}=x_{50},\\ \unicode[STIX]{x1D6FC}_{7}=x_{70},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

in which we denote by $x_{ij}$ the constant coefficients in $\unicode[STIX]{x1D6FC}_{i}$ in front of $\hat{I} _{j}$ (with the exception of $x_{i0}$ , which comes alone). We recall that $\hat{I} _{5}=\text{Tr}[\hat{\unicode[STIX]{x1D640}}^{2}]$ , $\hat{I} _{7}=\text{Tr}[\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}}]$ and $\hat{I} _{9}=\text{Tr}[\unicode[STIX]{x1D64C}\boldsymbol{\cdot }\hat{\unicode[STIX]{x1D640}}^{2}]$ .

As previously indicated, in simple shear flow, not all of the $x_{ij}$ can be separately tested against numerics. In the non-redundant componentwise basis (4.2), the linear model corresponds to functions $P_{-}$ , $P_{+}$ and $P_{12}$ of the form

(5.3) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}P_{-}=a_{-},\\ P_{+}=a_{+}+b_{+}Q_{+}+c_{+}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},\\ P_{12}=a_{12}+b_{12}Q_{+}+c_{12}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where we have introduced seven independent constants (the $a$ , $b$ and $c$ ).

Interestingly, this model is strongly constrained for the dynamics of $Q_{-}$ . More precisely, we now establish that equations (4.2) and (5.3) imply a rigorous condition on the post-reversal behaviour of $Q_{-}$ . After reversal, $\unicode[STIX]{x1D64C}$ reaches a steady-state value $\unicode[STIX]{x1D64C}^{+SS}$ whose components, from (4.2) and (5.3), satisfy

(5.4) $$\begin{eqnarray}\displaystyle 0=a_{-}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{-}^{+SS}+2Q_{12}^{+SS}. & & \displaystyle\end{eqnarray}$$

Subtracting this from the first line of (4.2) and defining $\unicode[STIX]{x0394}\unicode[STIX]{x1D64C}:=\unicode[STIX]{x1D64C}-\unicode[STIX]{x1D64C}^{+SS}$ , we get

(5.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\dot{Q}_{-}=a_{-}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})\unicode[STIX]{x0394}Q_{-}+2\unicode[STIX]{x0394}Q_{12}. & & \displaystyle\end{eqnarray}$$

Using the property of the pre-reversal steady state that $Q_{-}(\unicode[STIX]{x1D6FE}=0)=Q_{-}^{-SS}=Q_{-}^{+SS}$ , this can be integrated as

(5.6) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}Q_{-}(\unicode[STIX]{x1D6FE})=2\int _{0}^{\unicode[STIX]{x1D6FE}}\text{e}^{a_{-}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})[\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FE}^{\prime }]}\unicode[STIX]{x0394}Q_{12}(\unicode[STIX]{x1D6FE}^{\prime })\,\text{d}\unicode[STIX]{x1D6FE}^{\prime }. & & \displaystyle\end{eqnarray}$$

From this, we can see in particular that, until $\unicode[STIX]{x0394}Q_{12}(\unicode[STIX]{x1D6FE})$ changes sign, the sign of $\unicode[STIX]{x0394}Q_{-}(\unicode[STIX]{x1D6FE})$ has to stay the same as $\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})\operatorname{sgn}(\unicode[STIX]{x0394}Q_{12}(0))$ , irrespective of the values of the $a$ , $b$ and $c$ coefficients. This constraint is a strong prediction for any linear model that can be easily compared with the numerical data.

Figure 5.

Simulation data for $\unicode[STIX]{x0394}Q_{-}$ and $\unicode[STIX]{x0394}Q_{12}$ against post-reversal strain for three volume fractions, (a) $\unicode[STIX]{x1D719}=0.4$ , (b) $\unicode[STIX]{x1D719}=0.5$ and (c) $\unicode[STIX]{x1D719}=0.55$ . Any fabric evolution model that is linear in the fabric itself predicts that $\unicode[STIX]{x0394}Q_{-}$ cannot change sign before $\unicode[STIX]{x0394}Q_{12}$ does, from (5.6). This condition is unambiguously violated for all of our simulation data. Insets: zoom-in on the region where the condition is violated.

Figure 5 shows the post-reversal evolution of $\unicode[STIX]{x0394}Q_{-}$ and $\unicode[STIX]{x0394}Q_{12}$ for the three volume fractions we investigated. In every case, $\unicode[STIX]{x0394}Q_{-}$ changes sign long before $\unicode[STIX]{x0394}Q_{12}$ . We moreover verified that this behaviour is observed for any cut-off length $\unicode[STIX]{x1D716}$ . This allows us to conclude that no linear model can yield a good fit to the data. This is the first important conclusion of our work. We are not aware of any similar rejection of linear fabric evolution models in previous work that does not address strain reversal.

Beyond ruling out the possibility of finding a suitable linear model, the simulations allow us to conduct diagnostics, identifying the source of the difficulty. Between $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FE}+\text{d}\unicode[STIX]{x1D6FE}$ , $\unicode[STIX]{x1D64C}$ evolves through three distinct processes: the advection of pairs of particles with separation $h<\unicode[STIX]{x1D716}$ ; the birth of near-contacts counted in $\unicode[STIX]{x1D64C}$ when a pair of particles has separation $h>\unicode[STIX]{x1D716}$ at strain $\unicode[STIX]{x1D6FE}$ but $h<\unicode[STIX]{x1D716}$ an interval $\unicode[STIX]{x1D6FE}+\text{d}\unicode[STIX]{x1D6FE}$ later; and the death of near-contacts in the opposite case. Writing $r^{b}$ and $r^{d}$ for respectively the rate of near-contact births and deaths per near-contact, $\unicode[STIX]{x1D64C}^{b}$ and $\unicode[STIX]{x1D64C}^{d}$ for the fabric of near-contact births and deaths, and $\dot{\unicode[STIX]{x1D64C}}^{a}$ for the rate of change of $\unicode[STIX]{x1D64C}$ due to the advection of near-contacts not instantaneously being born or dying, we have (see appendix C)

(5.7) $$\begin{eqnarray}\displaystyle \dot{\unicode[STIX]{x1D64C}}=\dot{\unicode[STIX]{x1D64C}}^{a}+r^{b}\unicode[STIX]{x1D64C}^{b}-r^{d}\unicode[STIX]{x1D64C}^{d}. & & \displaystyle\end{eqnarray}$$

Instead of considering directly $\dot{\unicode[STIX]{x1D64C}}$ as a function of $\unicode[STIX]{x1D64C}$ , we can consider each term $\dot{\unicode[STIX]{x1D64C}}^{a}$ , $r^{b}\unicode[STIX]{x1D64C}^{b}$ and $r^{d}\unicode[STIX]{x1D64C}^{d}$ separately as functions of $\unicode[STIX]{x1D64C}$ . Trying to fit, using linear least squares, $\dot{\unicode[STIX]{x1D64C}}^{a}$ as a function of $\unicode[STIX]{x1D64C}$ from simulation data with the linear model in (4.2), we obtain an excellent result, shown in figure 6(a–c). This implies that the need for a nonlinear model is not due to the advective part of the evolution, but rather due to the birth and death of near-contacts. This is confirmed in figure 6(d–i), showing that the linear model fails for both the birth and death contributions, and gives particularly bad results when the birth and death contributions take their largest values in amplitude, i.e. when good accuracy matters most.

Figure 6.

The three $\dot{\unicode[STIX]{x1D64C}}$ contributions, $\dot{\unicode[STIX]{x1D64C}}^{a}$ (a–c), $r^{b}\unicode[STIX]{x1D64C}^{b}$ (d–f) and $r^{d}\unicode[STIX]{x1D64C}^{d}$ (g–i), as defined by (5.7), as functions of strain. For each contribution the three independent components $\dot{Q}_{-}^{a}$ , $\dot{Q}_{+}^{a}$ and $\dot{Q}_{12}^{a}$ (respectively, $r^{b}\unicode[STIX]{x1D64C}_{-}^{b}$ , $r^{b}\unicode[STIX]{x1D64C}_{+}^{b}$ and $r^{b}\unicode[STIX]{x1D64C}_{12}^{b}$ , and $r^{d}\unicode[STIX]{x1D64C}_{-}^{d}$ , $r^{d}\unicode[STIX]{x1D64C}_{+}^{d}$ and $r^{d}\unicode[STIX]{x1D64C}_{12}^{d}$ ) are shown respectively in (a,d,g), (b,e,h) and (c,f,i) (black lines). Each component is compared to a fit to the linear model (5.3) (colour lines).

5.2 Quadratic Hand equation

Having ruled out the possibility of an adequate linear model, we consider the general quadratic model

(5.8) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}P_{-}=a_{-}+b_{-}Q_{+}+c_{-}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},\\ P_{+}=a_{+}+b_{+}Q_{+}+c_{+}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}+d_{+}Q_{+}^{2}+e_{+}Q_{12}^{2}+f_{+}Q_{-}^{2},\\ P_{12}=a_{12}+b_{12}Q_{+}+c_{12}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}+d_{12}Q_{+}^{2}+e_{12}Q_{12}^{2}+f_{12}Q_{-}^{2}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

The parameter space of this model, with 15 dimensions, is difficult to explore. The strategy we adopt can be found in § 5.1 of Cheng, Jia & Wang (Reference Cheng, Jia and Wang2007), and is to numerically differentiate $\unicode[STIX]{x1D64C}$ (see appendix C), and use linear least squares to fit these quadratic models, with $\unicode[STIX]{x1D64C}$ taken from simulation data, to $\dot{\unicode[STIX]{x1D64C}}$ . The logic behind this is that, if the model is capable of yielding a good qualitative fit to $\unicode[STIX]{x1D64C}$ , it will also provide an adequate fit to $\dot{\unicode[STIX]{x1D64C}}$ , particularly if the fit is quantitatively good as well.

This approach very often fails for quadratic models governed by (5.8), because the parameter set obtained in this way may, when used to evolve $\unicode[STIX]{x1D64C}$ from its initial condition, cause $\unicode[STIX]{x1D64C}$ to grow unbounded. It is easy to see why: quadratic models correspond to the overdamped dynamics of a three-dimensional vector $\boldsymbol{s}=\{Q_{-},Q_{+},Q_{12}\}$ in a cubic potential ${\dot{s}}_{i}=-\unicode[STIX]{x1D735}_{s_{i}}[A_{jkl}s_{j}s_{k}s_{l}+B_{jk}s_{j}s_{k}+C_{k}s_{k}]$ , with $\unicode[STIX]{x1D63C}$ , $\unicode[STIX]{x1D63D}$ and $\unicode[STIX]{x1D63E}$ tensors depending on the coefficients of (5.8). A cubic potential is generically not confining, so that unless the initial conditions are within the basin of attraction of a local potential minimum, $\boldsymbol{s}$ will grow unbounded. Finding a well-behaved quadratic fabric evolution model then involves identifying parameters for which $\{Q_{-}^{-SS},Q_{+}^{-SS},Q_{12}^{-SS}\}$ lies within such a basin. This appears to be an unsatisfactory approach at best.

By setting some parameters of (5.8) to zero, it is possible to restrict the parameter space enough to find a model with a bounded evolution. For instance, the best quadratic model we could find takes the form

(5.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}P_{-}=a_{-}+b_{-}Q_{+}+c_{-}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12},\\ P_{+}=a_{+}+b_{+}Q_{+}+c_{+}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}+e_{+}Q_{12}^{2},\\ P_{12}=a_{12}+b_{12}Q_{+}+c_{12}\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}+e_{12}Q_{12}^{2}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

This contains 11 free parameters, which we fit initially by linear least squares on (4.2), that is, minimizing over the parameters $a_{-}$ , $b_{-}$ , $c_{-}$ , $a_{+}$ , $b_{+}$ , $c_{+}$ , $e_{+}$ , $a_{12}$ , $b_{12}$ , $c_{12}$ and $e_{12}$ the quantity

(5.10) $$\begin{eqnarray}\displaystyle X_{1} & = & \displaystyle \int _{\unicode[STIX]{x1D6FE}=0}^{\unicode[STIX]{x1D6FE}=3}\{(\dot{Q}_{-}^{data}-P_{-}[Q_{+}^{data},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}^{data},Q_{-}^{2\,data}]\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{-}^{data}-2Q_{12}^{data})^{2}\nonumber\\ \displaystyle & & \displaystyle +\,(\dot{Q}_{+}^{data}-P_{+}[Q_{+}^{data},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}^{data},Q_{-}^{2\,data}]\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}}))^{2}\nonumber\\ \displaystyle & & \displaystyle +(\dot{Q}_{12}^{data}-P_{12}[Q_{+}^{data},\operatorname{sgn}(\dot{\unicode[STIX]{x1D6FE}})Q_{12}^{data},Q_{-}^{2\,data}]+{\textstyle \frac{1}{2}}Q_{-}^{data})^{2}\!\}\,\text{d}\unicode[STIX]{x1D6FE}\end{eqnarray}$$

with strain series $\unicode[STIX]{x1D64C}^{data}$ and $\dot{\unicode[STIX]{x1D64C}}^{data}$ taken from the simulation data. Note that there is some freedom in choosing the upper $\unicode[STIX]{x1D6FE}$ limit in the integral of (5.10); a larger maximum $\unicode[STIX]{x1D6FE}$ will favour fits that are accurate in steady state, while a smaller one will favour good accuracy in the early phase. The value $\unicode[STIX]{x1D6FE}=3$ picked here is the one that we found to give the most balanced results between the two extremes. After this initial fit, we further apply a gradient descent to minimize the linear least squares between the actual model prediction $\unicode[STIX]{x1D64C}^{fit}$ (obtained by integrating equation (4.2) over strain with initial conditions from the simulation data $\unicode[STIX]{x1D64C}^{data}(\unicode[STIX]{x1D6FE}=0)$ ) and the numerical data $\unicode[STIX]{x1D64C}^{data}$ , that is, minimizing

(5.11) $$\begin{eqnarray}\displaystyle X_{2}=\int _{\unicode[STIX]{x1D6FE}=0}^{\unicode[STIX]{x1D6FE}=3}[(Q_{-}^{fit}-Q_{-}^{data})^{2}+(Q_{+}^{fit}-Q_{+}^{data})^{2}+(Q_{12}^{fit}-Q_{12}^{data})^{2}]\,\text{d}\unicode[STIX]{x1D6FE} & & \displaystyle\end{eqnarray}$$

to get final best-fit parameters.

Depending on the component of interest, this procedure leads to fits ranging from good to excellent, as shown in figure 7. The only obvious weakness is a moderate but increasing discrepancy for $Q_{-}$ when $\unicode[STIX]{x1D719}$ increases. This is probably due to the smaller $r^{b}$ and $r^{d}$ at low volume fraction (see figure 8), yielding flatter, easier-to-fit curves.

Figure 7.

Fits of the quadratic model described by (5.9) to simulation data for the three components of $\unicode[STIX]{x1D64C}$ (respectively, $Q_{-}$ , $Q_{+}$ and $Q_{12}$ from a–i) against post-reversal strain for $\unicode[STIX]{x1D719}=0.4$ (a,d,g), $\unicode[STIX]{x1D719}=0.5$ (b,e,h) and $\unicode[STIX]{x1D719}=0.55$ (c,f,i).

Figure 8.

Simulation data for the fractional birth rate $r^{b}$ (blue) and death rate $r^{d}$ (red) against post-reversal strain for (a) $\unicode[STIX]{x1D719}=0.4$ , (b) $\unicode[STIX]{x1D719}=0.5$ and (c)  $\unicode[STIX]{x1D719}=0.55$ .

Because of the difficulty to find a well-behaved model, we cannot perform a thorough exploration of parameter space, and there may be a better-fitting quadratic model than (5.9). In any case the fitting power of quadratic models is, in three dimensions, in principle satisfactory when it comes to matching the numerically observed fabric evolution. However, we believe this to be largely coincidental, rather than faithful to the underlying physics, and stemming simply from the large number of parameters available. This is partly because there is no reason to expect the parameters of a mechanistically faithful model to need fine-tuning to avoid the generic blow-up described previously. However, a stronger reason to reject such models is found by making a detour to the two-dimensional case, as we describe next.