A general constitutive model for dense, fine particle suspensions validated in many geometries

TL;DR

This paper introduces a comprehensive three-dimensional continuum model for dense, fine particle suspensions that captures complex shear thickening, jamming fronts, and non-Newtonian behaviors across various flow geometries.

Contribution

It develops a novel, predictive, and general mixture theory-based model with a micro-structural variable, validated against multiple flow scenarios and behaviors.

Findings

Accurately captures shear thickening behaviors (CST and DST).

Models propagation of shear and impact jamming fronts.

Represents non-Newtonian effects like 'running on oobleck'.

Abstract

Fine particle suspensions (such as cornstarch mixed with water) exhibit dramatic changes in viscosity when sheared, producing fascinating behaviors that captivate children and rheologists alike. Recent examination of these mixtures in simple flow geometries suggests inter-granular repulsion is central to this effect --- for mixtures at rest or shearing slowly, repulsion prevents frictional contacts from forming between particles, whereas, when sheared more forcefully, granular stresses overcome the repulsion allowing particles to interact frictionally and form microscopic structures that resist flow. Previous constitutive studies of these mixtures have focused on particular cases, typically limited to two-dimensional, steady, simple shearing flows. In this work, we introduce a predictive and general, three-dimensional continuum model for this material, using mixture theory to couple the fluid and particle phases. Playing a central role in the model, we introduce a micro-structural state variable, whose evolution is deduced from small-scale physical arguments and checked with existing data. Our space- and time-dependent model is implemented numerically in a variety of unsteady, non-uniform flow configurations where it is shown to accurately capture a variety of key behaviors: (i) the continuous shear thickening (CST) and discontinuous shear thickening (DST) behavior observed in steady flows, (ii) the time-dependent propagation of `shear jamming fronts', (iii) the time-dependent propagation of `impact activated jamming fronts', and (iv) the non-Newtonian, `running on oobleck' effect wherein fast locomotors stay afloat while slow ones sink.

Tables

Table 1. Table 1: Material parameters for curves in figure 1 . μ c = 1 subscript 𝜇 𝑐 1 \mu_{c}=1 for all fits. μ h subscript 𝜇 ℎ \mu_{h} determined from a 𝑎 a and ϕ m subscript italic-ϕ 𝑚 \phi_{m} (see [ 4 ] ).

Fit	${\varphi }_c$	${\varphi }_j$	$a_0$	$a_{\infty }$	${\tau }^{\ast }/\left({\eta }_0{\dot{\gamma }}_0\right)$	${\eta }_0/{\eta }_B$
(a)	0.585	0.65	0.455	0.841	1.74	0
(b)	0.585	0.65	0.455	0.841	1.44	0.254$\left({\varphi }_j-\varphi \right)$

Table 2. Table 2: Material parameters for curves in figure 2 . μ c = 0.3 subscript 𝜇 𝑐 0.3 \mu_{c}=0.3 for [ 12 ] , [ 17 ] and [ 20 ] . μ c = 1.19 subscript 𝜇 𝑐 1.19 \mu_{c}=1.19 for [ 16 ] and [ 22 ] . Δ = 0 Δ 0 \Delta=0 for all experiments.

Fit	${\varphi }_c$	${\varphi }_j$	$a_0$	$a_{\infty }$	${\tau }^{\ast }$ [Pa]	${\eta }_0$ [mPa$\cdot$s]
[12]	0.57	0.65	1.08	1.08	5.21	56
[17]	0.548	0.642	1.08	1.08	17.93	130
[20](a)	0.56	0.617	1.577	1.577	423.0	2.4
[20](b)	0.563	0.610	1.09	1.538	3.5	2.4
[16]	0.38	0.454	1.00	1.00	0.116	0.89
[22]	0.436	0.535	0.8	0.8	5.35	6

Table 3. Table 3: Parameters of η ^ B ​ ( ϕ ) subscript ^ 𝜂 𝐵 italic-ϕ \hat{\eta}_{B}(\phi) for curves in figure 2 .

Fit	${\phi }_1$ [(Pa$\cdot$s)-1]	${\alpha }_1$	${\phi }_2$ [(Pa$\cdot$s)-1]	${\alpha }_2$
[12]	$8.83\times {10}^6$	6	0	0
[17]	$1.79\times {10}^{11}$	12	$4.5$	0.1
[20](a)	$3\times {10}^6$	5	$9\times {10}^{-2}$	0.1
[20](b)	$1.0\times {10}^6$	6	0	0
[16]	$4.0\times {10}^4$	3	0	0
[22]	$1.86\times {10}^4$	4	$9.3\times {10}^{-3}$	0.1

Table 4. Table 4: Material parameters for simulations shown in figure 3 .

${\varphi }_0$	${\eta }_0$ [mPa$\cdot$s]	$\rho$ [kg/m3]	$G$ [Pa]	$K$ [Pa]	$K_0$
0.505	10	1620	$3.8\times {10}^5$	$8.3\times {10}^5$	0.06

Table 5. Table 5: Material parameters for simulations shown in figures 4 and 5 .

${\eta }_0$ [mPa$\cdot$s]	${\rho }_s$ [kg/m3]	${\rho }_f$ [kg/m3]	$G$ [Pa]	$K$ [Pa]	$K_0$
7 or 1	1620	1000	$3.8\times {10}^7$	$8.3\times {10}^7$	1.0

Table 6. Table 6: Material parameters for simulations shown in figures 6 and 7 .

${\eta }_0$ [mPa$\cdot$s]	${\rho }_s$ [kg/m3]	${\rho }_f$ [kg/m3]	$G$ [Pa]	$K$ [Pa]	$K_0$
4.6	1620	1620	$3.8\times {10}^7$	$8.3\times {10}^7$	0.2

Table 7. Table 7: Material parameters for simulations shown in figures 8 and 9 .

${\eta }_0$ [mPa$\cdot$s]	${\rho }_s$ [kg/m3]	${\rho }_f$ [kg/m3]	$G$ [Pa]	$K$ [Pa]	$K_0$
1.0	1620	1620	$3.8\times {10}^5$	$8.3\times {10}^5$	1.0

Table 8. Table 8: Summary of equations for full three-dimensional model.

Rule	Expression
Granular Phase Mass Conservation	$\frac{D^s{\overline{\rho }}_s}{Dt}+{\overline{\rho }}_s\sum _{i=1}^3\frac{\partial v_s^{}{}_i{}^{}}{\partial x_i}=0$
Fluid Phase Mass Conservation	$\frac{D^f{\overline{\rho }}_f}{Dt}+{\overline{\rho }}_f\sum _{i=1}^3\frac{\partial v_f^{}{}_i{}^{}}{\partial x_i}=0$
Fluid Phase True Density	$\frac{n}{{\rho }_f}\frac{D^f{\rho }_f}{Dt}=-\sum _{i=1}^3\frac{\partial \left(\left(1-n\right)v_s^{}{}_i{}^{}+nv_f^{}{}_i{}^{}\right)}{\partial x_i}$
Granular Phase Momentum Balance	${\overline{\rho }}_s\frac{D^{s}v_s^{}{}_i{}^{}}{Dt}={\overline{\rho }}_s{g}_i-f_d^{}{}_i{}^{}-\left(1-n\right)\frac{\partial p_f}{\partial x_i}+\sum _{i=1}^3\frac{\partial {\tilde{\sigma }}_{ij}}{\partial x_j}$
Fluid Phase Momentum Balance	${\overline{\rho }}_f\frac{D^{f}v_f^{}{}_i{}^{}}{Dt}={\overline{\rho }}_f{g}_i+f_d^{}{}_i{}^{}-n\frac{\partial p_f}{\partial x_i}+\sum _{i=1}^3\frac{\partial \tau _f^{}{}_{ij}{}^{}}{\partial x_j}$
Darcy’s Drag Law	$f_d^{}{}_i{}^{}=\frac{18\varphi \left(1-\varphi \right){\eta }_0}{d^2}\hat{F}(\varphi ,\mathrm{Re})\left(v_s^{}{}_i{}^{}-v_f^{}{}_i{}^{}\right)$
Fluid Phase Pore Pressure	$p_f=\kappa \ln \left(\frac{{\rho }_f}{{\rho }_{0f}}\right)$
Fluid Phase Shear Stress	$\tau _f^{}{}_{ij}{}^{}={\eta }_0\left(1+\frac{5}{2}\varphi \right)\left(\frac{\partial v_f^{}{}_i{}^{}}{\partial x_j}+\frac{\partial v_f^{}{}_j{}^{}}{\partial x_i}-\frac{2}{3}\sum _{k=1}^3\frac{\partial v_f^{}{}_k{}^{}}{\partial x_k}{\delta }_{ij}\right)$
Granular Phase Effective Stress	${\dot{\tilde{\sigma }}}_{ij}=2G\left(D_{ij}-{\tilde{D}}_{ij}^p\right)+\left(K-\frac{2}{3}G\right)\left(D_{kk}-{\tilde{D}}_{kk}^p\right){\delta }_{ij}+W_{ik}{\tilde{\sigma }}_{kj}-{\tilde{\sigma }}_{ik}{W}_{kj}$
Granular Phase Strain-Rate and Spin Tensors	$D_{ij}=\frac{1}{2}\left(\frac{\partial v_s^{}{}_i{}^{}}{\partial x_j}+\frac{\partial v_s^{}{}_j{}^{}}{\partial x_i}\right),W_{ij}=\frac{1}{2}\left(\frac{\partial v_s^{}{}_i{}^{}}{\partial x_j}-\frac{\partial v_s^{}{}_j{}^{}}{\partial x_i}\right)$
Granular Phase Plastic Flow Rate	${\tilde{D}}_{ij}^p=\frac{{\dot{\overline{\gamma }}}^p}{2\overline{\tau }}\left({\tilde{\sigma }}_{ij}+\tilde{p}{\delta }_{ij}\right)+\frac{1}{3}\left(\beta {\dot{\overline{\gamma }}}^p+{\dot{\xi }}_1+{\dot{\xi }}_2\right){\delta }_{ij}$
Equivalent Granular Shear Stress	$\overline{\tau }=\sqrt{\frac{1}{2}\left({\sigma }_{kl}+\tilde{p}{\delta }_{kl}\right)\left({\sigma }_{kl}+\tilde{p}{\delta }_{kl}\right)}$
Granular Pressure	$\tilde{p}=-\frac{1}{3}{\tilde{\sigma }}_{kk}$
Dilation Angle	$\beta =K_3\left(\varphi -{\varphi }_{eq}\right)$
Critical State Packing Fraction	${\varphi }_{eq}=\frac{{\varphi }_m}{1+a{I}_m}$
Frictional model for ${\varphi }_m$	${\varphi }_m={\varphi }_j+\left({\varphi }_c-{\varphi }_j\right)f$
Frictional model for $a$	$a=a_0+\left(a_{\infty }-a_0\right)f$
Evolution of Granular Structure	$\dot{f}=K_0{\dot{\overline{\gamma }}}^p{\left(\frac{\overline{\tau }}{{\tau }^{\ast }}\right)}^{3/2}\left(f_m-f\right)-K_0\left({\dot{\overline{\gamma }}}^p+\frac{\overline{\tau }}{{\eta }_B}+{\dot{\xi }}_{ϵ}\right)f$
Internal Friction Coefficient	${\mu }_p={\mu }_1+\frac{{\mu }_2-{\mu }_1}{1+\left(b/I_m\right)}+\frac{5}{2}\left(\frac{\varphi I_v}{a{I}_m}\right)$
Inertial Numbers	$I_v=\frac{{\eta }_0{\dot{\overline{\gamma }}}^p}{\tilde{p}},I={\dot{\overline{\gamma }}}^{p}d\sqrt{\frac{{\rho }_s}{p}},I_m=\sqrt{I^2+2I_v}$
Granular Shear Flow Rule
Granular Separation Rule
Granular Compaction Rule

Equations

$$\frac{\tau}{\eta_{0}\dot{\gamma}}=\eta_{r}=1+\frac{5}{2}\phi\bigg{(}\frac{\phi_{m}}{\phi_{m}-\phi}\bigg{)}+2\mu_{c}\bigg{(}\frac{a\phi}{\phi_{m}-\phi}\bigg{)}^{2},$$

$$\frac{\tilde{p}}{\eta_{0}\dot{\gamma}}=\eta_{n}=2\bigg{(}\frac{a\phi}{\phi_{m}-\phi}\bigg{)}^{2}$$

$${\varphi }_m={\hat{\varphi }}_m(f)={\varphi }_j+({\varphi }_c-{\varphi }_j)f,$$

$$a=\hat{a}(f)=a_0+(a_{\infty }-a_0)f,$$

$$\frac{\dot{f}}{K_0\dot{\gamma }}=H(f_m-f)-Sf,$$

$$H=\bigg{(}\frac{\bar{\tau}}{\tau^{*}}\bigg{)}^{3/2},$$

$$S=1+\frac{\bar{\tau }}{{\eta }_B\dot{\gamma }}+\frac{{\dot{\xi }}_{ϵ}}{\dot{\gamma }},$$

$${\hat{f}}_m(\varphi )=\{\begin{array}{ll}1,& \text{if }\varphi \le {\varphi }_c\\ \frac{\varphi -{\varphi }_j}{{\varphi }_c-{\varphi }_j},& \text{if }{\varphi }_c<\varphi \le {\varphi }^{\ast }\\ \frac{{\varphi }^{\ast }-{\varphi }_j}{{\varphi }_c-{\varphi }_j},& \text{if }\varphi >{\varphi }^{\ast },\end{array}$$

$$\eta_{B}=\hat{\eta}_{B}(\phi)=\bigg{(}\sum_{i=1}^{i_{\text{max}}}\varphi_{i}(\phi_{j}-\phi)^{\alpha_{i}}\bigg{)}^{-1}$$

$${\tilde{D}}_{ij}^p=\frac{{\dot{\bar{\gamma }}}^p}{2\bar{\tau }}({\tilde{\sigma }}_{ij}+\tilde{p}{\delta }_{ij})+\frac{1}{3}(\beta {\dot{\bar{\gamma }}}^p+{\dot{\xi }}_1+{\dot{\xi }}_2){\delta }_{ij},$$

$$\bar{\tau }=\sqrt{\frac{1}{2}({\tilde{\sigma }}_{kl}+\tilde{p}{\delta }_{kl})({\tilde{\sigma }}_{kl}+\tilde{p}{\delta }_{kl})},\tilde{p}=-\frac{1}{3}{\tilde{\sigma }}_{kk},$$

$$\frac{\dot{f}}{K_0\dot{\gamma }}=H(f_m-f)-Sf$$

$$\dot{f}=0⟹f=\frac{H{f}_m}{H+S}.$$

$$H=\hat{H}(\bar{\tau }/{\tau }^{\ast }),$$

$$S=\bigg{(}1+\frac{\dot{\xi}_{B}}{\dot{\gamma}}+\frac{\dot{\xi}_{\epsilon}}{\dot{\gamma}}\bigg{)},$$

$${\dot{\xi }}_B={\dot{\xi }}_{ϵ}=0,\dot{f}=0⟹f=f(\bar{\tau })=\frac{\hat{H}(\bar{\tau }/{\tau }^{\ast })}{\hat{H}(\bar{\tau }/{\tau }^{\ast })+1}.$$

$$H=\hat{H}(\bar{\tau}/\tau^{*})=\bigg{(}\frac{\bar{\tau}}{\tau^{*}}\bigg{)}^{3/2}.$$

$${\varphi }_W=\frac{m_{cs}}{m_{cs}+m_l},$$

$$\varphi =\frac{(1+\lambda )(1-\beta )(m_{cs}/{\rho }_{cs})}{(1-\beta )(m_{cs}/{\rho }_{cs})+(m_l/{\rho }_l)+\beta (m_{cs}/{\rho }_w)},$$

$$\varphi =\frac{(1+\lambda )(1-\beta )}{(1-\beta )+(\frac{1-{\varphi }_W}{{\varphi }_W})({\rho }_{cs}/{\rho }_l)+\beta ({\rho }_{cs}/{\rho }_w)},$$

$$\frac{\tau}{\eta_{0}\dot{\gamma}}=\eta_{r}=1+\frac{5}{2}\phi\bigg{(}\frac{\phi_{m}}{\phi_{m}-\phi}\bigg{)}+2\mu_{c}\bigg{(}\frac{a\phi}{\phi_{m}-\phi}\bigg{)}^{2},$$

$${\varphi }_m={\hat{\varphi }}_m(f)={\varphi }_j+({\varphi }_c-{\varphi }_j)f,$$

$$a=\hat{a}(f)=a_0+(a_{\infty }-a_0)f.$$

$$f\to 0\quad\implies\quad\eta_{r}(\phi)=1+\frac{5}{2}\phi\bigg{(}\frac{\phi_{j}}{\phi_{j}-\phi}\bigg{)}+2\mu_{c}\bigg{(}\frac{a_{0}\phi}{\phi_{j}-\phi}\bigg{)}^{2}.$$

$$f\to 1\quad\implies\quad\eta_{r}(\phi)=1+\frac{5}{2}\phi\bigg{(}\frac{\phi_{c}}{\phi_{c}-\phi}\bigg{)}+2\mu_{c}\bigg{(}\frac{a_{\infty}\phi}{\phi_{c}-\phi}\bigg{)}^{2}.$$

$${\hat{f}}_m(\varphi )=\frac{{\varphi }_j-\varphi }{{\varphi }_j-{\varphi }_c},\text{if }{\varphi }_c\le \varphi \le {\varphi }^{\ast }.$$

$$\dot{f}=0⟹f=\frac{H}{H+S}\frac{{\varphi }_j-\varphi }{{\varphi }_j-{\varphi }_c},\text{if }{\varphi }_c\le \varphi \le {\varphi }^{\ast },$$

$${\varphi }_m={\varphi }_j+(\varphi -{\varphi }_j)\frac{H}{H+S},\text{if }{\varphi }_c\le \varphi \le {\varphi }^{\ast }.$$

$$\frac{\bar{\tau}}{\eta_{0}\dot{\gamma}}=\frac{5}{2}\bigg{(}\frac{\phi^{2}}{\phi_{j}-\phi}\bigg{)}\bigg{(}\frac{H+S}{S}\bigg{)}+2\mu_{c}\bigg{(}\frac{a\phi}{\phi_{j}-\phi}\bigg{)}^{2}\bigg{(}\frac{H+S}{S}\bigg{)}^{2},\text{\quad if \quad}\phi_{c}\leq\phi\leq\phi^{*}.$$

$$\lim_{\bar{\tau}\to\infty}{H}=\bigg{(}\frac{\bar{\tau}}{\tau^{*}}\bigg{)}^{3/2},$$

Full text

1 Introduction

The behavior of granular materials suspended in fluid media has been a major topic of study for over a century. These types of mixtures are present in many industrial, geotechnical, and biological engineering problems, spanning length scales from tens of meters to millimeters. Of particular interest in this work is the behavior of chemically stable, hard, frictional particles suspended in viscous fluids as found in industrial processes and studied in soil mechanics.

Though thoroughly studied and classified, a unifying constitutive model relating mixture stresses, strains, and strain rates across all material types and flow regimes remains elusive (see [38]). One challenge for such a model is capturing the non-Newtonian, shear thickening behavior that is often observed in mixtures like water-cornstarch and water-poly(methyl methacrylate) when the mean particle diameter, $d$, is less than $\sim$10 $\mu$m. At low volume fractions, the apparent viscosity of such mixtures grows steadily with increasing shearing rate (CST); however, at high volume fractions, the apparent viscosity of these mixtures can jump several orders of magnitude with very little change in measured shearing rate (DST).

Recent theoretical work in [42], experimental observations in [20], [11], [15], and [8], and simulations reported in [35], [27], [28], and [36] have shown that the shear thickening behavior observed in these mixtures is a direct result of inter-granular repulsion and its effect on the dilation behavior of these mixtures. In the presence of relatively small applied stresses, the particles in the mixture will interact through lubrication forces in the suspending medium and behave like a granular material with low internal friction. In the presence of relatively large applied stresses, the particles in the mixture will be forced into frictional contact, drastically increasing the apparent internal friction coefficient. In dense suspensions, this frictional transition can cause the granular skeleton to dilate.

Prior work modeling general fluid-sediment mixtures in [4] has focused on hard, frictional, non-Brownian, non-repulsive particle suspensions ($d\gtrsim 100$ $\mu$m for common engineering mixtures). The modeling framework proposed in that work combines the empirical relations presented in [38], [7] and [1] with the two-phase mixture theories developed in [13] and [23] into a single constitutive theory. In this work, we build upon this framework to produce a more general model for granular suspensions that accounts for the effect of inter-granular repulsion on mixtures with mean grain diameters $d\lesssim 100$ $\mu$m. Using a custom numerical scheme, we validate this model in transient and inhomogeneous flows, in both two-dimensional and three-dimensional geometries. The model is shown to replicate the unusual landmark features of these shear-thickening suspensions such as the ability to run across them but not walk.

2 Model

General Theory

In continuum modeling of fluid-particle mixtures, we consider the materials that constitute the mixture independently. The individual solid particles (or grains) are homogenized into a single continuum body called the granular phase which has a density $\bar{\rho}_{s}=\phi\rho_{s}$ and velocity components ${v_{s}}_{i}$. Here $\phi$ defines the volume fraction of the granular phase and $\rho_{s}$ defines the density of each individual grain. Similarly, the fluid which fills the space between the grains is homogenized into a single continuum body called the fluid phase with density $\bar{\rho}_{f}=(1-\phi)\rho_{f}$ and velocity components ${v_{f}}_{i}$. Here $\rho_{f}$ is called the true density of the fluid. A thorough theoretical description of this process of homogenization can be found in [4] and [23].

From the basic laws of mass conservation, momentum and energy balance, and entropy imbalance, it is possible to define a complete set of governing equations for the evolution of density, velocity, and stress within each phase. These equations account for (i) momentum exchange between the phases through a buoyant force and a Darcy-like inter-phase drag, (ii) dilation and contraction of the granular phase through a specialized elastic-plastic constitutive model, and (iii) separation of mixture stresses into components from the granular skeleton, $\tilde{\sigma}_{ij}$, and components from the pore fluid, ${\tau_{f}}_{ij}$ and $p_{f}$. These equations are combined to form the model proposed in [4]. This model is shown to be accurate in a wide range of flow regimes and geometries; however, it does not account for the effects of inter-granular repulsion and corresponding structure evolution. In this work, we show that reformulating two previously constant parameters from the model shown in [4], $a$ and $\phi_{m}$, is sufficient to capture the physics of particle repulsion and accurately model steady CST and DST as well as several transient and dynamic behaviors observed in the dense, fine grain suspensions we seek to model.

The full set of equations that define our proposed three-dimensional, time-dependent, two-phase model can be found in the Supplemental Information; however, several key features of this model are easily illustrated by considering its behavior in steady two-dimensional (or quasi-two-dimensional) shearing flows. This behavior is described by the relationship between the steady mixture shearing rate $\dot{\gamma}$, the mixture shear stress $\tau$, and the granular phase pressure $\tilde{p}$ . For flows dominated by viscous effects (e.g. water-cornstarch; $(\rho_{s}d^{2}\dot{\gamma})/\eta_{0}\ll 1$), that steady state relationship is identical to the model proposed in [31] and presented in [7] and is expressed as follows,

[TABLE]

[TABLE]

with $\rho_{s}$ the density of the grains, $\eta_{0}$ the fluid phase viscosity, $\phi$ the granular volume fraction, and $\mu_{c}$ the flow resistance due to frictional granular contacts. $\phi_{m}$ and $a$ are dilation parameters and have a significant effect on the mixture stress and granular phase pressure.

The parameter $\phi_{m}$ limits the range of volume fractions where steady flow is possible (i.e. for $\phi\geq\phi_{m}$, shearing the material results in continuous growth of the shear stress, sometimes referred to as a jammed state) and the parameter $a$ modifies the critical state volume fraction (see $\phi_{eq}$ in [4] and Supplemental Information) and associated Reynolds’ dilation. $\phi_{m}$ and $a$ are often treated as material constants; however, as proposed in [42] and fully defined in [36], $\phi_{m}$ and $a$ are better described as functions with $\phi_{m}$ bounded between $\phi_{j}$ and $\phi_{c}$ (with $\phi_{j}>\phi_{c}$) as follows,

[TABLE]

and $a$ bounded between $a_{0}$ and $a_{\infty}$ as follows,

[TABLE]

with $f$ a scalar measure of the anisotropy within the granular structure and the fraction of frictional granular interactions (bounded between 0 and 1). Conceptually, the material parameters $\phi_{j}$ and $\phi_{c}$ correspond to the range of volume fractions where steady shearing flow can occur. For $\phi<\phi_{c}$, steady flow can occur at all shearing rates. For $\phi\geq\phi_{j}$, no steady flow is possible; the mixture is jammed. For $\phi_{c}\leq\phi<\phi_{j}$, steady flow is only possible at some shearing rates. $a_{0}$ and $a_{\infty}$ are fitting parameters.

Unlike in previous models where $f$ was a function of granular stress directly, we propose that $f$ evolves over time according to the following rule,

[TABLE]

with $K_{0}$ a dimensionless parameter, $H$ a dimensionless hardening rate, $S$ a dimensionless softening rate, and $f_{m}$ an upper bound on the steady value of $f$ given by $f_{m}=\hat{f}_{m}(\phi).$ Such an upper bound on $f$ was first proposed for low volume fractions in [34] and is an important feature in this model, especially as it relates to rheo-chaos (see [22] and [19]).

The first term in Eq. (5), $H\ (f_{m}-f)$, describes the evolution of $f$ toward $f_{m}$ and, by Eq. (1) and Eq. (3), the associated increase in effective viscosity, $\eta_{r}$. Based on the work of [27] (see Supplemental Information), we propose the following form for $H$,

[TABLE]

with $\tau^{*}$ a repulsive stress scale with units of Pa (see $\tau_{\text{min}}$ in [9], and $\sigma^{*}$ in [36]), $\bar{\tau}=(\mu_{c}+\mu_{h})\tilde{p}$ the effective granular shear stress, and $\mu_{h}$ the flow resistance due to fluid-mediated granular interactions (see [7]). This form of $H$ captures the strengthening of the granular phase as the granular stress forcing grains together overcomes inter-granular repulsion. The second term in Eq. (5), $-Sf$, describes the evolution of $f$ toward 0 and the associated decrease in effective viscosity, $\eta_{r}$. We propose the following form for $S$,

[TABLE]

with $\eta_{B}$ a measure of the resistance of the granular skeleton to buckling induced degradation with units of Pa$\cdot$s, and $\dot{\xi}_{\epsilon}$ a rate scale for Brownian or electrostatic repulsion with units s*-1*. The three terms in $S$ reflect three proposed mechanisms which will break down granular structure. First, macroscopic shearing of the mixture will cause grains to slip past each other, altering their structure (see [10] for discussion of ‘plastic rearrangement’ in jammed materials).

Second, we postulate that force chains can undergo load-initiated reorganization ($\bar{\tau}/\eta_{B}\dot{\gamma}$) presumably through chain buckling, inducing structural degradation without macroscopic deformation. In systems of grains with large but finite elastic stiffness $G_{p}$, the small relative contact area $A_{c}/d^{2}$ of the particles means the buckling load $\tau_{c}\sim G_{p}A_{c}^{2}/d^{4}$ need not be large. The proposed form might then be seen as the product of a statistical measure of critical chains ($\sim\bar{\tau}/\tau_{c}$) and a time-scale of buckling degradation given by the properties of the mixture ($t_{c}\tau_{c}/\eta_{0}=\hat{\varphi}(\rho G_{p}d^{2}/\eta_{0}^{2},A_{c}/d^{2},\phi)$, i.e. $\eta_{B}\propto\tau_{c}t_{c}$). And lastly, Brownian or electrostatically driven diffusion will drive grains apart ($\dot{\xi}_{\epsilon}/\dot{\gamma}$); in flowing materials at room temperature, this contribution to $\dot{f}$ should be negligible.

Steady Behavior of Rheologically Stable Flows

We begin analyzing the proposed model by considering the simulated three-dimensional shearing flows in [36]. These flows were composed of between 500 and 2000 stiff, spherical particles undergoing simple shear in a periodic domain. The discrete particle interactions were modeled using standard contact laws in conjunction with additional forces to account for lubrication and inter-granular repulsion. The steady state relationship between $\dot{\gamma}/\dot{\gamma}_{0}$ and $\tau/(\eta_{0}\dot{\gamma}_{0})$ (with $\eta_{0}\dot{\gamma_{0}}\propto\tau^{*}$) reported in that work (see figure 1) contains several regimes where $\partial\dot{\gamma}/\partial\tau\leq 0$. These regimes can lead to unstable behavior under imposed shearing rates (see [29]) and imposed shear stresses (see [19]). The latter instability manifests in rheologically chaotic behavior (i.e. transient inhomogeneities in the flow occur; rheo-chaos). The lack of reported rheo-chaos in [36] suggests that the behavior of larger systems of particles may differ at packing fractions near and above $\phi_{c}$.

We can, however, still examine the steady state response of the model in Eq. (5) for such rheologically stable shearing flows. To do this, we first define a functional form for $\hat{f}_{m}(\phi)$. The results in [27] are generated using the method shown in [36] and suggest that $f_{m}$ is insensitive to changes in $\phi$. We therefore let $f_{m}=1$ for such mixtures. With $f_{m}$ determined, we calculate and fit the steady shear response of our model to the data in [36] by solving for $\dot{f}=0$ at different volume fractions $\phi$ and shear rates $\dot{\gamma}/\dot{\gamma}_{0}$ (with $\tau/\eta_{0}\dot{\gamma}_{0}$ determined from Eq. (1)) as shown in figure 1. The relevant material parameters for these fits are provided in table 1. (See Supplemental Information for specific details about our model fitting procedure.)

Pseudo-Steady Behavior of Rheologically Chaotic Flows

We continue analyzing the proposed model by considering experimental analysis of repulsive grain mixtures. Although many authors have explored the phenomenology of these mixtures, complete characterization of their behavior is complicated by several factors which require additional modeling considerations beyond the ideal case discussed in the prior section. The scaling of $\tau^{*}$ with respect to mean grain diameter $d$ is often reported as $\tau^{*}\propto d^{-2}\text{ or }d^{-3}$ (see [17], [20], and [9]) and is observed to have vanishing effect on the response of common granular suspensions with $d\gtrsim 100$ $\mu$m; this puts practical limitations on the minimum system size (quantified by number of grains) that can be analyzed. Slip along boundaries (see [16] and [33]) and mixture breakdown are observed when applied shear stresses exceed approximately $10^{4}$ Pa, limiting the range of mixture responses that can be probed. Additionally, starches and other porous particles swell significantly in suspension, obfuscating the true volume fraction of the mixture (see [22]).

Nevertheless, experimentation is the most direct way to examine the behavior of these particle suspensions and has yielded many important observations. Among these is the rheo-chaos reported in [22], [12], and [6]; these rheologically chaotic flows are characterized by large, rapid changes in the measured mixture shearing rate at constant applied stresses and measured relationships between $\dot{\gamma}$ and $\bar{\tau}$ that generally have $\partial\dot{\gamma}/\partial\bar{\tau}\geq 0$. More thorough study is necessary to understand the mechanics and transient effects of rheo-chaos, but if we consider the long-time average behavior of the shearing flows reported in the literature to be representative of the average local behavior, we can examine the steady state response of the model in Eq. (5) for such rheologically chaotic flows.

We begin this analysis by defining a new expression for $\hat{f}_{m}(\phi)$ based on the experimental observations reported in [22] and [16] as follows,

[TABLE]

with $\phi^{*}=\phi_{j}+(\phi_{c}-\phi_{j})\Delta$ and $\Delta$ an implicit function of the compliance of the granular skeleton. By observation, we expect $\Delta\to 0$ for large systems of irregular grains (where rheo-chaos is observed) and $\Delta\to 1$ for small systems of stiff, spherical grains (where rheo-chaos is not observed). In the latter case, $f_{m}\to 1$ as in the idealized mixtures of the previous section.

Given this form of $\hat{f}_{m}(\phi)$ it is possible to determine a critical shearing rate associated with DST, $\dot{\gamma}_{\text{DST}}$, which is the limiting value of $\dot{\gamma}$ as $\bar{\tau}$ increases for the volume fractions in the range $\phi_{c}\leq\phi\leq\phi^{*}$ (see Supplemental Information). Based on the distribution of apparent $\dot{\gamma}_{\text{DST}}$ values in the literature and the trends of $\dot{\gamma}_{c}$ reported in [3], we propose the following form of the structural resistance to buckling induced degradation $\eta_{B}$,

[TABLE]

with $\varphi_{i}$ representing material parameters with units of (Pa$\cdot$s)-1 and $\alpha_{i}$ representing dimensionless scaling factors.

With these expressions determined, we calculate and fit the steady shear response of our model to experimental results reported in [12], [17], [20], [16], and [22] by solving for $\dot{f}=0$ while satisfying the relationship in Eq. (1) at at the relevant volume fractions and shearing rates as shown in figure 2. All data from the various experiments are seen to be well represented by the generic form proposed. The parameters for these steady shearing experiments can be found in tables 2 and 3. The wide range of $\tau^{*}$ values reported in table 2 is likely linked to the scaling of $\tau^{*}\propto d^{-2}$ or $d^{-3}$ and changes in surface chemistry between different granular materials (see the effects of hydrogen bonding on ‘shear jamming’ and DST in [24]). [12], [17], and [20] performed experiments using poly(methyl methacrylate) beads with diameters $d$ between roughly 0.3$\mu$m and 4$\mu$m while [16] and [22] performed experiments using cornstarch with reported grain diameters often between 5$\mu$m and 20$\mu$m. (See Supplemental Information for specific details about fitting $\dot{\gamma}_{\text{DST}}$ and determining $\phi$ from [22].)

Complete, Time-Dependent Expression of Model

The steady state measurement of DST and CST is only one part of the puzzle of repulsive grain suspensions. As noted in [40], a fully coupled constitutive model is necessary to capture the observed dynamic behavior of these mixtures as epitomized by the ‘running on oobleck’ effect. This effect describes the ability of a person to run over the surface of a mixture of standard cornstarch and water without sinking (within a certain range of volume fractions); however, if a person walks over the same mixture, they rapidly sink to the bottom.

In the next section, we will show that the model proposed in this work is capable of capturing this behavior; however, we must first express the time-dependent, three-dimensional form of our model. The complete set of governing equations (many of which are unchanged from the model developed in [4]) can be found in the Supplemental Information. Here, we focus only on the granular phase stress components $\tilde{\sigma}_{ij}$ which obey a Maxwell-like model. We assume that the deformation rate tensor, $D_{ij}=\frac{1}{2}(\partial v_{si}/\partial x_{j}+\partial v_{sj}/\partial x_{i})$, can be decomposed into an elastic part $\tilde{D}^{e}_{ij}$ and plastic part $\tilde{D}^{p}_{ij}$ with ${v_{s}}_{i}$ the components of the granular phase velocity. This allows us to express the rate of change of the granular phase stress $\tilde{\sigma}_{ij}$ as a function of the elastic deformation rate, the granular shear modulus $G$ and bulk modulus $K$. The components of the plastic part of the deformation rate tensor can be expressed using the form given in [4] as follows,

[TABLE]

with,

[TABLE]

$\delta_{ij}$ the Kronecker delta function, and $\dot{\bar{\gamma}}^{p}$ the scalar (equivalent) plastic shear rate that generalizes $\dot{\gamma}$ to unsteady, three-dimensional flows. The expression in Eq. (10) captures the two mechanisms of granular flow: shearing and dilation. The second term in Eq. (10) (representing dilation of the granular phase) captures the effects of Reynolds dilation ($\beta\dot{\bar{\gamma}}^{p}$), pure expansion ($\dot{\xi}_{1}$), and pure compaction ($-\dot{\xi}_{2}$). $\beta$, $\dot{\bar{\gamma}}^{p}$, and $\dot{\xi}_{2}$, are functions of the granular stress $\tilde{\sigma}_{ij}$, $a$, and $\phi_{m}$; $\dot{\xi}_{1}$ is a function of inter-granular cohesion only. Implementation of the changes made in our proposed model is achieved by replacing $a$ and $\phi_{m}$ in the functions for $\beta$, $\dot{\bar{\gamma}}^{p}$, and $\dot{\xi}_{2}$ with Eq. (3) and Eq. (4), and substituting $\dot{\bar{\gamma}}^{p}$ in place of $\dot{\gamma}$ in Eq. (5).

3 Results

In this section, we show that the fully three-dimensional, time-dependent, two-phase form of our model accurately reproduces the observed transient behavior of cornstarch-water mixtures in annular shear and impact experiments as well as the ‘running on oobleck’ effect. To perform this analysis, we implement our model in the numerical framework shown in [4]. This numerical framework is an adaptation of the material point method (MPM) for simulating mixtures as two overlapping continua. This version of MPM is very similar to the method shown in [2] and differs from the original method described in [39] in that each continuum phase of the mixture is represented by independent sets of material point tracers. These material point tracers act as quadrature points for solving the weak form of the momentum balance equations on a static background simulation grid. In addition, these tracers move with the continua that they represent. To understand the results shown in this work, it suffices to understand that the velocities and strain-rates of the continua are represented on the nodes of the simulation grid while displacements and stresses are represented on the material point tracers.

As in [4], we assume that the granular phase of the mixture is elastically stiff (i.e. elastic moduli $G$ and $K\gg\bar{\tau}$) so that any significant deformation of the material comes from plastic flow ($\dot{\bar{\gamma}}^{p}$). In all of the simulations reported in this work, we choose values of the fluid bulk modulus and the moduli $G$ and $K$ to be large enough such that any effects of elasticity do not alter the flow solution (to do this, we require that all elastic deformations are smaller than 1%) but not so large that they significantly impact the quality of the numerical calculations. In this way, we ensure that our results are representative of elastically stiff mixtures while also avoiding well known numerical artifacts in MPM (see kinematic locking in [30] and the ringing instability in [18]) that affect the smoothness of the spatial stress field. Further discussion of this point can be found in the Supplemental Information along with a brief note about the validity of the elastically stiff assumption.

Dynamic Shear Jamming

Here we show that our model can capture the reported propagation of ‘shear jamming fronts’ in annular shear experiments; in particular, we are interested in the starkly different mixture responses reported at different applied shearing rates. One example of this phenomena can be found in [33] which reports several experiments involving density matched mixtures of cornstarch, water, gylcerol, and CsCl. In the annular shear cell described in that work, an applied inner wall velocity of $u_{i}=0.84$ mm/s produces shearing flow that rapidly approaches the diffusive flow profile expected from a Newtonian fluid; however, when the applied inner wall velocity is multiplied by a factor of 100, the mixture begins ‘jamming’ at the inner wall. This ‘jammed’ region is characterized by nearly rigid rotation (see figure 3) and a well defined boundary that propagates outward with a near constant rate.

We recreate the reported experiment with our continuum model using an axisymmetric form of MPM on a 36$\times$50 element Cartesian grid with 1mm$\times$1mm resolution. Our model is calibrated to the experiments using the steady flow parameters for the cornstarch-water/glycerol experiments of [22] (see tables 2 and 3) and material parameters as in table 4. The only parameters that were tuned to fit experimental observations were $K_{0}$ and $\phi_{0}$. The values chosen for $\eta_{0}$, $d$, and $\rho$ are informed by the materials used in the mixtures ($d$ in all simulations is chosen to be 5.85 $\mu$m; see [37]). $K_{3}$ and $K_{4}$ are chosen according to the findings of [4] (4.715 and 0, respectively). The simulation results in figure 3 show excellent agreement with experimental observations neglecting wall slip effects that were not implemented in our model.

Impact Activated Solidification

In this section, we show that our model accurately captures the solidification of constarch-water/glycerol mixtures under impact. Similar to the ‘shear jamming’ observed in the previous section, this ‘impact activated solidification’ results in two strikingly different material responses. At high impact speeds the mixture responds like a solid, sometimes causing the impactor to bounce (see [40]); at low impact speeds the mixture responds like a fluid, flowing around the intruder with ease. In addition, as predicted in [40] and directly observed in [32] and [21], at high impact speeds, a ‘solid plug’ forms at the impact site and propagates outward (another type of ‘jamming front’).

We begin analyzing the response of our (fully dynamic) model under impact by recreating several of the experiments reported in [40] where an aluminum rod is launched at a mixture of cornstarch, water, and glycerol. Measurement of the rod motion was taken visually and with an accelerometer, and mixture displacement fields were visualized using X-ray imaging. We calibrate our model to the experiments using the steady flow parameters for [22] in tables 2 and 3 and material parameters in table 5 ($d$, $K_{3}$, and $K_{4}$ are chosen as in the previous section). The simulations are performed using axisymmetric MPM on Cartesian grids with 1mm$\times$1mm resolution. We only tune $K_{0}$ to accurately reproduce the displacement fields (see figure 4) and rod dynamics (see figure 5) reported in [40]. By utilizing the full continuum momentum balance, our model naturally captures observed added mass effects from material inertia.

The behavior of our model during impact shows remarkable similarity to experiments. The stiff response of the mixture at high impact speeds (as shown by the acceleration curves in figure 5) is contrasted by the steady sinking that is observed once the rods have slowed (as shown by the velocity curves in figure 5). The X-ray imaging performed in [40] allowed observation of the final displacement field within the mixture (see figure 4) and the results are indicative of a ‘solidification front’ propagating during impact.

We continue the analysis of our model during impact by recreating the controlled intrusion experiments described in [21] with a focus on measuring the rate of propagation of the ‘solidification front’ as the intruder enters the mixture. In these experiments, a rod was driven with a constant velocity $v_{0}$ into a density matched cornstarch-water/glycerol/CsCl mixture, and the internal flow of the mixture was visualized using high speed ultrasound imaging. We calibrate our model to the experiments using the parameters for [22] found in tables 2 and 3 and material parameters in table 6 ($d$, $K_{3}$, and $K_{4}$ are chosen as before). The simulations are run using axisymmetric MPM on a 50$\times$35 element Cartesian grid with 1mm$\times$1mm resolution.

Tuning $K_{0}$, we are able to closely match the reported flow fields (see figure 6) from [21]. As the intruder enters the mixture, a ‘solid plug’ of material forms below the intruder and grows outward in all directions. This ‘solid plug’ is characterized by nearly rigid vertical motion with a well defined boundary that grows with near constant rate. The rate of propagation of this boundary in the transverse $v_{ft}$ and longitudinal $v_{fl}$ directions is often quite different (as shown in figure 6(f) and 6(g)); however, these rates appear to obey a particular scaling rule that is independent of volume fraction: $(v_{fl}-v_{0})/v_{ft}\approx 2$. Remarkably, our simple three-dimensional model naturally captures this reported scaling (see figure 7; results are shown in terms of the transverse and longitudinal front propagation factors $k_{t}=v_{ft}/v_{0}$ and $k_{l}=v_{fl}/v_{0}-1$, respectively). This close agreement of model predictions to experimental observations is unmatched by any other models in the literature.

Running on Oobleck

In this section we demonstrate that our model can capture the ‘running on oobleck’ effect by simulating a spoked elastic wheel driving over the surface of a density matched cornstarch-fluid mixture similar to those described above. The simulation is run using two-dimensional MPM on a 150$\times$75 element (2m$\times$1m) periodic Cartesian grid using the steady flow parameters for [22] in tables 2 and 3 with material parameters in table 7. ($d$, $K_{3}$, and $K_{4}$ are chosen as in the previous sections). The four spoked elastic wheel used in this simulation has radius 15cm and density 3000 kg/m3 and is driven with constant rotation rate $\omega$ starting on the mixture’s surface (see figures 8 and 9).

As shown in figures 9 and 8 respectively, the wheel driving at $\omega=1.0$ s*-1* will rapidly sink while the wheel driving at $\omega=10.0$ s*-1* will run along the surface. This difference in behavior is caused by the same mechanism explored in the previous section. In the slow case ($\omega=1.0$ s*-1*), the spokes of the wheel are traveling too slowly to produce a solid-like response; the entire wheel slips into the mixture as it would in a fluid. However, in the fast case ($\omega=10.0$ s*-1*), when the spokes of the wheel impact the mixture surface, they cause rapid solidification of the surrounding mixture — see the spacial $f$ field in figures 8 and 9 (plotted using a post-processing routine from [14]). This allows that spoke to support the weight of the wheel until the next spoke strikes the surface further on; the wheel runs along the top of the mixture.

4 Conclusion

The model proposed in this work shows remarkable accuracy in predicting the behavior of repulsive particle suspensions in both steady and unsteady flows. By combining a model for the evolution of granular micro-structure with the governing equations for fluid-sediment mixtures presented in [4], we have created a single constitutive model for chemically stable, hard, frictional particles suspended in viscous fluids. In addition, implementation of this model in the numerical framework proposed in [4] allows simulation of the coupled behavior of particles and fluid through a wide range of particle sizes, volume fractions, and flow regimes (including independent motion of each phase and interaction with solid bodies). We have demonstrated that this model captures the landmark features of these mixtures including (i) DST and CST, (ii) propagation of ‘shear jamming fronts’, (iii) the propagation of ‘impact activated jamming fronts’, and (iv) the ‘running on oobleck’ effect.

There are a number of areas remaining for further research and improvement of this model. In this work, we have focused on the time-average behavior of rheologically chaotic flows; however, the basic model may be used to predict the time-accurate behavior of these flows through a fabric tensor contribution to the evolution of $f$ (as noted in [12] and [38]). In addition, future work will need to address the dependence of $\tau^{*}$, $\eta_{B}$, and $K_{0}$ on changes in particle and fluid material properties such as $\eta_{0}$, $d$, $\rho$, and $\theta$ (with $\theta$ temperature; see [41]) and the additional dependence of $K_{0}$ on stress scale $\bar{\tau}$ and shearing rate $\dot{\gamma}$. In the Supplemental Information, we describe a simple relaxation experiment which could help determine these dependencies. Numerically, we have shown that our model can be implemented in MPM (see [4]) and accurately reproduce the time-dependent behavior of these highly non-linear mixtures in non-trivial geometries. Although our initial tests are promising, further work should also be done to curb known issues in the numerical framework such as kinematic locking (see [30]) and the ringing instability (see [18]), which can cause spurious artifacts in certain parameter ranges if not dealt with carefully.

5 Supplemental Information

Complete Table of Equations

The complete, three-dimensional, two-phase model proposed in this work is an extension of the model presented in [4]. As described in the main document, this model changes the definition of the dilation parameters $a$ and $\phi_{m}$ from material constants (as in [4]) to functions of the granular micro-structure $f$. A summary of the equations that define the complete model can be found in table 8.

These equations govern the behavior of the two independent continua that form the mixture: the granular phase and the fluid phase. The mass conservation rules govern the evolution of the fluid phase true density $\rho_{f}$ and effective density $\bar{\rho}_{f}=n\rho_{f}$ as well as the granular phase effective density $\bar{\rho}_{s}=\phi\rho_{s}$. Here we assume that $\rho_{s}$ is constant (incompressible grains) and define $n$, the mixture porosity, as equal to $(1-\phi)$. The two momentum balance equations describe the evolution of the components of the fluid and granular phase velocities ${v_{f}}_{i}$ and ${v_{s}}_{i}$ according to gradients in the fluid phase pore pore pressure $p_{f}$, components of the effective granular stress tensor $\tilde{\sigma}_{ij}$, and components of the fluid phase shear stress tensor ${\tau_{f}}_{ij}$. In addition, these rules also define the components of the inter-phase Darcy’s Law drag that acts between the two material phases, ${f_{d}}_{i}$; the scalar function $\hat{F}(\phi,\textrm{Re})$ is defined in [4] and taken from [5].

The fluid phase stresses are modeled using a viscous, barotropic constitutive model which allows for weak compressibility and accounts for Einstein’s effective viscosity model for dilute suspensions.

The granular phase effective stress is modeled as an elastic-plastic solid. Through careful definition of the components of the plastic flow rate tensor $\tilde{D}^{p}_{ij}$ (given by the scalar rates $\dot{\bar{\gamma}}^{p}$, $\dot{\xi}_{1}$ and $\dot{\xi}_{2}$), this model can accurately describe stresses in the granular skeleton of a mixture at rest, in steady flow, and all regimes in between. In addition, the form of the plastic flow rate tensor in table 8 captures the different modes of deformation which a granular material can undergo: pure shear ($\dot{\bar{\gamma}}^{p}$), shear dilation ($\beta\dot{\bar{\gamma}}^{p}$), pure dilation ($\dot{\xi}_{1}$), and pure compaction ($-\dot{\xi}_{2}$). The scalar flow rates $\dot{\bar{\gamma}}^{p}$, $\dot{\xi}_{1}$ and $\dot{\xi}_{2}$ are themselves functions of the granular phase effective stress and are defined implicitly by the yield functions $f_{1}$, $f_{2}$, and $f_{3}$ respectively. These yield functions are written such that the granular phase captures the $\mu(I)$ dry inertial rheology from [25], the $\mu(I_{v})$ viscous inertial rheology from [7], and the $\mu(I_{m},I_{v})$ mixed inertial rheology from [1].

Determining $H$ from [27]

The proposed model for the evolution of granular micro-structure $f$ in dense suspensions has the following form,

[TABLE]

with H a hardening parameter, S a softening parameter, $f_{m}$ a limiting value for $f$ (taken here to be 1), and $K_{0}$ a scaling term. This expression for $\dot{f}$ has a zero which corresponds to steady behavior of the granular mixture,

[TABLE]

In this work, we have posited that

[TABLE]

with $\tau^{*}$ a repulsive stress scale and $\bar{\tau}$ a measure of the applied granular stress such that the rate of hardening of the granular micro-structure depends on the magnitude of the applied granular stress overcoming inter-granular repulsion. Additionally, we have posited that the softening behavior can be described as follows,

[TABLE]

with a constant term for the rate of structural decay in shear, $\dot{\xi}_{B}$ the rate of structural breakdown due to buckling of force chains, and $\dot{\xi}_{\epsilon}$ the rate of structural breakdown due to diffusion.

If $\dot{\xi}_{B}$ and $\dot{\xi}_{\epsilon}$ are taken to be negligible in comparison to the shearing rate $\dot{\gamma}$, then we find that the steady behavior of $f$ is a function of stress only,

[TABLE]

In dense suspensions, the granular stress $\bar{\tau}$ will account for most of the total mixture stress $\tau$ such that $f(\bar{\tau})\approx f(\tau)$. This direct dependence of $f$ on the mixture stress $\tau$ is in agreement with the models proposed in [42] ($f(\tau)\approx 1-e^{-\tau/\tau^{*}}$) and [36] ($f(\tau)\approx e^{-\tau^{*}/\tau}$) as well as the simulated flows reported in [27]. As shown in figure 10, there is strong agreement between these models and our prediction for $f(\bar{\tau})$ when

[TABLE]

Estimating $\phi$ from $\phi_{W}$ in [22]

Due to the significant swelling of starch particles suspended in fluid, the mixture mass fractions $\phi_{W}$ is reported for the suspensions studied in [22] instead of the true volume fractions $\phi$ with,

[TABLE]

and $m_{cs}$ the mass of cornstarch in the mixture and $m_{l}$ the mass of suspending fluid in the mixture. To convert this to solid volume fraction $\phi$ for model fitting, we use the method described in [33] that attempts to account for pore space and particle swelling as follows,

[TABLE]

with $\lambda\approx 0.3$, $\beta\approx 0.1$, $\rho_{cs}$ the density of a cornstarch grain, $\rho_{l}$ the density of the suspending fluid, and $\rho_{w}$ the density of water. This expression can be written equivalently as follows,

[TABLE]

and is used to determine the values of $\phi$ shown in this work.

Fitting $a_{0}$ and $\phi_{j}$ to Steady Shearing Data

Here we describe a method for determining $a_{0}$ and $\phi_{j}$ from steady shearing data. Recall the expression for the effective viscosity of a mixture in steady shearing flow from [4],

[TABLE]

and the expressions for $a$ and $\phi_{m}$ from [36],

[TABLE]

[TABLE]

If we consider the behavior of a mixture at relatively low shearing rates and stresses, we expect that $f\to 0$. In this limit, we find an expression for effective viscosity $\eta_{r}$ that is identical to that from [7] with $\phi_{m}=\phi_{j}$ and $a=a_{0}$,

[TABLE]

We can fit this expression to the low viscosity measurements pulled from experiments (see figure 11) to determine the best values of $a_{0}$ and $\phi_{j}$. (Note that $\mu_{c}$ is chosen to reflect the dry granular behavior of the granular particles and is not fit using these curves.)

Fitting $a_{\infty}$ and $\phi_{c}$ to Steady Shearing Data

As in the previous section, we can fit $a_{\infty}$ and $\phi_{j}$ to experimental data if we consider the behavior of a mixture at relatively high shearing rates and stresses where we expect that $f\to 1$. In this limit, we find an expression for effective viscosity $\eta_{r}$ that is identical to that from [7] with $\phi_{m}=\phi_{c}$ and $a=a_{\infty}$,

[TABLE]

We can fit this expression to the high viscosity measurements pulled from experiments (see figure 12) to determine the best values of $a_{\infty}$ and $\phi_{c}$. (Note that we limit the range of volume fractions considered to flow curves which have a well defined secondary plateau.)

Fitting $\dot{\gamma}_{\text{DST}}$ to Steady Shearing Data

For rheologically stable flows, the parameters $\tau^{*}$ and $\eta_{B}$ are fit by observation with $\eta_{B}$ generally taken to be very large (see fits for [36] in main document). For rheologically chaotic flows, large applied shear stresses at packing fractions above $\phi_{c}$ appear to cause unbounded DST (mixture breakdown and wall slip limit the range of stresses that can be experimentally measured). To fit $\eta_{B}$ to such flows, we begin by recalling the proposed form of $\hat{f}_{m}(\phi)$ for volume fractions in the range of $\phi_{c}\leq\phi\leq\phi^{*}$,

[TABLE]

Combining this expression with Eq. (13), we find,

[TABLE]

and together with the expression for $\phi_{m}$ from [36],

[TABLE]

Substituting this expression into the an expression for the effective granular viscosity in [4], we have,

[TABLE]

We continue by recalling the forms of $H$ and $S$ proposed in this work; from these, it can be shown that,

[TABLE]

and,

[TABLE]

All of this combines together to give the following result,

[TABLE]

which can be reduced to an equation of $\dot{\gamma}$ and $\phi$ only to yield the expression for $\dot{\gamma}_{\text{DST}}$ shown in the main document,

[TABLE]

We can then fit this expression to the experimentally measured shearing rates associated with observed DST to determine a reasonable form for $\eta_{B}$ (see figure 13). Observation of the trends of $\dot{\gamma}_{\text{DST}}$ in the literature (and $\dot{\gamma}_{c}$ in [3]) suggest the form proposed for $\eta_{B}$ in the main document,

[TABLE]

such that the resistance to buckling induced degradation, $\eta_{B}$, grows as the volume fraction $\phi$ increases toward $\phi_{j}$ at a rate controlled by $\varphi_{i}$ and $\alpha_{i}$.

Ringing Instability in Cornstarch-Water Simulations using MPM

The material point method (MPM; see [39]) is a numerical scheme for solving dynamic mechanics problems and excels when history dependent materials undergo large deformations. Unlike standard Lagrangian methods, a static background grid is used to represent the material velocities and weak form test functions; using a fixed grid avoids the problems of mesh distortion that often accompany large material deformations. Unlike standard Eulerian methods, a fixed set of material point tracers are used to represent history dependent properties (e.g. $\tilde{\sigma}_{ij}$, $f$, $\phi$) and act as quadrature points for weak form integration; tracking history dependent quantities at fixed points in the material avoids the errors associated with numerical advection schemes. These two features of MPM lead to an unfortunate side effect which is common to many ‘particle methods’: the ringing instability (see [18]).

The ringing instability is an accumulated error in the material state represented on the material point tracers. This error is allowed to accumulate in the null space of the material-point-to-grid mapping matrix. At each time-step, the material state as represented on the material point tracers must be ‘mapped’ or integrated in order to update the grid representation of velocity. In order to maintain high accuracy when integrating the weak form equations, we often require that the number of points $N_{p}$ is greater than the number of grid coefficients $N_{n}$ (often equal to the number of nodes); this guarantees the existence of a null space in which error can grow. In order to minimize the effect of this error on the quality of the results presented in this work, we have generally chosen elastic moduli $G$ and $K$ (and the fluid bulk modulus $\kappa$) which are ‘large enough’ to avoid significant elastic deformations, but not ‘so large’ as to produce significantly spurious stress fields (see figure 14 where example pressure fields, $\tilde{p}$, beneath the elastic wheel in the ‘running on oobleck’ simulations are shown for three different sets of elastic moduli).

A Method for Determining $G$ and $K_{0}$

In the model presented in this work, we have assumed that the material which constitutes the granular phase of the mixture is elastically stiff. This assumption follows from two observations: (i) the individual grains in mixtures which exhibit DST have elastic moduli on the order of GPa (e.g. poly(methyl methacrylate)) and (ii) the largest stresses measured in these mixtures are often smaller than a few MPa. If the elastic moduli $G$ and $K$ of the bulk granular material are of a similar order to the elastic moduli of the individual grains, then it is reasonable to expect that $G$ and $K\gg\bar{\tau}$. For this reason, we have neglected to find realistic values of $G$ and $K$ for use in our model and have assumed that most of the observed dynamics of these mixtures are dominated by plastic flow ($\dot{\bar{\gamma}}^{p}$) instead of elastic deformation.

In this section, we describe the behavior of our proposed model as it relates to stress relaxation in the mixture after the cessation of simple, quasi-two-dimensional shearing flow for the purposes of experimentally determining the values of $G$ and $K_{0}$. For reference, the basic equation of our model can be found in table 8. If the granular material is sheared to steady-state, then it can be shown that $\phi<\phi_{m}$, $\dot{f}=0$, and,

[TABLE]

Given a steady shear stress measurement of $\bar{\tau_{0}}$ (and assuming $\phi_{c}$, $\phi_{j}$, $a_{0}$, $a_{\infty}$, $\tau^{*}$, and $\eta_{B}$ have already been determined), it is then possible to find the steady value of $f$, $f_{0}$.

With the steady response of the system known, we can then consider its behavior when the material is stopped suddenly; that is, we are interested in the time-accurate measurement of the shear stress, $\bar{\tau}(t)$, after the applied shearing rate is set to zero, $D_{ij}=W_{ij}=0$. In this regime, we define $\dot{\bar{\gamma}}^{p}$ according to Eq. (35), $\dot{f}$ according to the rules in table 8, and $\dot{\bar{\tau}}$ as follows (derived from expressions in table 8):

[TABLE]

We can then integrate Eq. (36) and the expression for $\dot{f}$ in table 8 to find,

[TABLE]

Example stress relaxation curves for $\bar{\tau}(t)$ at varying values of $G$ and $K_{0}$ can be found in figure 15. Prior work on examining the relaxation characteristics of these mixtures can be found in [12] and [26]. An important note about the results found in those works is that the characteristic stress relaxation time is on the order of tens of seconds for some of the mixtures considered. In order for our model to reproduce those relaxation times (see figure 15(a)), the granular shear modulus $G$ needs to be on the order of 10 kPa, much smaller than what we would normally expect. This time-scale difference suggests that the elastic behavior of these mixtures is more complicated than simple linear elasticity and further work on our model will be necessary to capture this behavior. However, the qualitative similarity between the curves shown in figure 15 and those found in [26] is a promising indication that we are capturing the correct phenomena.