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

See allHide authors and affiliations

Edited by David A. Weitz, Harvard University, Cambridge, MA, and approved August 27, 2019 (received for review May 23, 2019)

## Significance

Shear-thickening suspensions display peculiar phenomena, exhibiting solid-like and fluid-like behaviors at the same volume fraction. These materials have puzzled researchers since the 1930s and have seen renewed interest due to advances in computational and experimental methods as well as new infrastructural and military applications that might exploit these behaviors. We introduce a continuum model for shear-thickening suspensions, which correctly predicts realistic 3D flows including transient and inhomogeneous cases. Our approach builds the microscopic physics of particle repulsion into a framework for general fluid–grain mixtures, bridging to the fields of geomechanics and granular rheology. Using a material point method routine, we provide a robust model implementation and demonstrate the method’s broader utility for a variety of fluid-particle flow problems.

## 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. Examination of these mixtures in simple flow geometries suggests intergranular repulsion and its influence on the frictional nature of granular contacts 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 2D, steady, simple shearing flows. In this work, we introduce a predictive and general, 3D 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 microstructural 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, nonuniform flow configurations where it is shown to accurately capture a variety of key behaviors: 1) the continuous shear-thickening (CST) and discontinuous shear-thickening (DST) behavior observed in steady flows, 2) the time-dependent propagation of “shear jamming fronts,” 3) the time-dependent propagation of “impact-activated jamming fronts,” and 4) the non-Newtonian, “running on oobleck” effect, wherein fast locomotors stay afloat while slow ones sink.

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 (e.g., transport of waste material, erosion of shorelines and riverbeds). 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.

Although thoroughly studied and classified, a unifying constitutive model relating mixture stresses, strains, and strain rates across all material types and flow regimes remains elusive (1). 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 ∼10 μ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).

Theoretical work in refs. 2 and 3, experimental observations in refs. 4⇓⇓⇓–8, and simulations reported in refs. 9⇓⇓–12 have shown that the shear-thickening behavior observed in these mixtures is a direct result of intergranular repulsion inhibiting the formation of frictional granular contacts and the resulting effect on the dilation behavior of the mixture. 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 interparticle friction. In the presence of relatively large applied stresses, the particles in the mixture will be forced into frictional contact, drastically increasing the bulk resistance to flow. In dense suspensions, this frictional transition can cause the granular skeleton to dilate. Further, recent observations reported in refs. 13⇓⇓–16 and simulations shown in ref. 17 have indicated that the growth and decay of these dilating structures in the granular skeleton plays an important role in the time-dependent behavior of these mixtures.

Prior work modeling general fluid–sediment mixtures in ref. 18 has focused on hard, frictional, non-Brownian, nonrepulsive particle suspensions (

## Model Development

### 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

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 1) momentum exchange between the phases through a buoyant force and a Darcy-like interphase drag, 2) dilation and contraction of the granular phase through a specialized elastic-plastic constitutive model, and 3) separation of mixture stresses into components from the granular skeleton, **1** and **2** below), 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 3D, time-dependent, 2-phase model can be found in *SI Appendix*, Table S1; however, several key features of this model are easily illustrated by considering its behavior in steady, constant volume, homogeneous, 2D (or quasi-2D) shearing flows. This behavior is described by the relationship between the steady mixture shearing rate

The parameter *SI Appendix*, Table S1) and associated Reynolds’ dilation. **1**) and causing “jamming” through a wider range of volume fractions. The material parameters

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

The first term in Eq. **5**, **1** and **3**, the associated increase in effective viscosity, *SI Appendix*, Fig. S1), we propose the following form for H,**5**, ^{−1}. The 3 terms in S reflect 3 proposed mechanisms that will break down granular structure and thereby decrease the fraction of grains in frictional contact. First, macroscopic shearing of the mixture will cause grains to slip past each other, altering their structure (see ref. 28 for discussion of “plastic rearrangement” in jammed materials). Second, we postulate that force chains can undergo load-initiated reorganization (

### Steady Behavior of Rheologically Stable Flows.

We begin analyzing the proposed model by considering the simulated 3D shearing flows in ref. 12. These flows were composed of between 500 and 2,000 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 intergranular repulsion. The steady-state relationship between

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 **1**), as shown in Fig. 1. The relevant material parameters for these fits are provided in Table 1. (See *SI Appendix* 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 that require additional modeling considerations beyond the ideal case discussed in the prior section. The scaling of *E*) is observed when applied shear stresses exceed ∼

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 refs. 25, 32, and 33; these rheologically chaotic flows are characterized by large, rapid changes in the measured mixture shearing rate at constant applied stresses and measured relationships between **5** for such rheologically chaotic flows.

We begin this analysis by defining an expression for

Given this form of *SI Appendix*). Based on the distribution of apparent ^{−1} and

With these expressions determined, we calculate and fit the steady-shear response of our model to experimental results reported in refs. 5, 25, and 30⇓–32 by solving for **1** at the relevant volume fractions and shearing rates as shown in Fig. 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 *SI Appendix* for specific details about fitting

### 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 ref. 36, 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, 3D form of our model. The complete set of governing equations (many of which are unchanged from the model developed in ref. 18) can be found in *SI Appendix*, Table S1. Here, we focus only on the granular-phase stress components **10** captures the 2 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 (**3** and **4** and substituting **5**.

## Results

In this section, we show that the fully 3D, time-dependent, 2-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 ref. 18. This numerical framework is an adaptation of the material point method (MPM) for simulating mixtures as 2 overlapping continua. This version of MPM is very similar to the method shown in ref. 37 and differs from the original method described in ref. 38 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 ref. 18, we assume that the granular phase of the mixture is “elastically stiff” (i.e., elastic moduli G and *SI Appendix*, 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 phenomenon can be found in ref. 13, which reports several experiments involving density-matched mixtures of cornstarch, water, glycerol, and CsCl. In the annular-shear cell described in that work, an applied inner wall velocity of

We recreate the reported experiment with our continuum model using an axisymmetric form of MPM on a 36 × 50 element Cartesian grid with 1 mm × 1 mm resolution. The upper surface of the simulated mixture is modeled as “free” (zero traction), the inner and outer surfaces are modeled using no-slip, rigid boundaries with assigned velocity (according to the relevant experimental parameters), and the bottom surface is modeled as a rigid wall with slip (zero friction). Our model is calibrated to the experiments using the steady-flow parameters for the cornstarch–water/glycerol experiments of ref. 25 (Tables 2 and 3) and material parameters as in Table 4. The only parameters that were tuned to fit experimental observations were

### Impact-Activated Solidification.

In this section, we show that our model accurately captures the solidification of cornstarch–water/glycerol mixtures under impact. Similar to the shear jamming observed in the previous section, this “impact-activated solidification” results in 2 strikingly different material responses. At high impact speeds, the mixture responds like a solid, sometimes causing the impactor to bounce (36); at low impact speeds, the mixture responds like a fluid, flowing around the intruder with ease. In addition, as predicted in ref. 36 and directly observed in refs. 14 and 42, 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 ref. 36 where an aluminum rod (

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 Fig. 5) is contrasted by the steady sinking that is observed once the rods have slowed (as shown by the velocity curves in Fig. 5). The X-ray imaging performed in ref. 36 allowed observation of the final displacement field within the mixture (Fig. 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 ref. 14 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

Tuning *F* and *G*); however, these rates appear to obey a particular scaling rule that is independent of volume fraction:

### 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 2D MPM on a 150 × 75 element (2 m × 1 m) periodic Cartesian grid using the steady-flow parameters for ref. 25 in Tables 2 and 3 with material parameters in Table 7 (d,

As shown in Figs. 8 and 9, respectively, the wheel driving at ^{−1} will rapidly sink, while the wheel driving at ^{−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 (^{−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 (^{−1}), when the spokes of the wheel impact the mixture surface, they cause rapid solidification of the surrounding mixture—see the spatial f field in Figs. 8 and 9 (plotted using a postprocessing routine from ref. 44). 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. Analysis of the stresses within the spoke of the fast wheel show peak vertical stresses of around 20 kPa after the spoke strikes the mixture, much larger than the approximate Archimedes pressure of 3 kPa due to the displacement of the mixture. (Video of these simulations can be found in Movies S1–S4.)

## 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 microstructure with the governing equations for fluid–sediment mixtures presented in ref. 18, 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 ref. 18 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 1) DST and CST, 2) propagation of shear jamming fronts, 3) the propagation of impact-activated jamming fronts, and 4) 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 refs. 1 and 32). In addition, future work will need to address the dependence of *SI Appendix*, we describe a simple relaxation experiment that could help determine these dependencies. Numerically, we have shown that our model can be implemented in MPM (18) and accurately reproduce the time-dependent behavior of these highly nonlinear mixtures in nontrivial 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 (39) and the ringing instability (40), which can cause spurious artifacts in certain parameter ranges if not dealt with carefully.

## Acknowledgments

We acknowledge support from NSF Grant CBET-1253228 and Army Research Office Grants W911NF-16-1-0440 and W911NF-15-1-0196.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: kkamrin{at}mit.edu.

Author contributions: A.S.B. and K.K. performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1908065116/-/DCSupplemental.

Published under the PNAS license.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- C. Clavaud,
- A. Bérut,
- B. Metzger,
- Y. Forterre

- ↵
- A. Fall,
- F. Bertrand,
- G. Ovarlez,
- D. Bonn

- ↵
- ↵
- ↵
- R. Mari,
- R. Seto,
- J. F. Morris,
- M. M. Denn

- ↵
- R. Mari,
- R. Seto,
- J. F. Morris,
- M. M. Denn

- ↵
- A. Singh,
- R. Mari,
- M. M. Denn,
- J. F. Morris

- ↵
- I. R. Peters,
- S. Majumdar,
- H. M. Jaeger

- ↵
- E. Han,
- I. R. Peters,
- H. M. Jaeger

- ↵
- R. Maharjan,
- E. Brown

- ↵
- R. Maharjan,
- S. Mukhopadhyay,
- B. Allen,
- T. Storz,
- E. Brown

- ↵
- O. Ozgen,
- M. Kallmann,
- E. Brown

- ↵
- A. S. Baumgarten,
- K. Kamrin

- ↵
- ↵
- L. Amarsid et al.

- ↵
- D. S. Drumheller

- ↵
- R. Jackson

- ↵
- J. F. Morris,
- F. Boulay

- ↵
- J. R. Royer,
- D. L. Blair,
- S. D. Hudson

- ↵
- M. Hermes et al.

- ↵
- M. Grob,
- A. Zippelius,
- C. Heussinger

- ↵
- ↵
- ↵
- R. Mari,
- R. Seto,
- J. F. Morris,
- M. M. Denn

- ↵
- W. J. Frith,
- P. d’Haene,
- R. Buscall,
- J. Mewis

- ↵
- A. Fall et al.

- ↵
- P. d’Haene,
- J. Mewis,
- G. G. Fuller

- ↵
- W. H. Boersma,
- P. J. M. Baets,
- J. Laven,
- H. N. Stein

- ↵
- ↵
- N. M. James,
- E. Han,
- R. Arturo Lopez de la Cruz,
- J. Jureller,
- H. M. Jaeger

- ↵
- ↵
- S. Bandara,
- K. Soga

- ↵
- ↵
- C. M. Mast,
- P. Mackenzie-Helnwein,
- P. Arduino,
- G. R. Miller,
- W. Shin

- ↵
- C. E. Gritton

- ↵
- M. Stasiak,
- M. Molenda,
- I. Opaliński,
- W. Błaszczak

- ↵
- I. R. Peters,
- H. M. Jaeger

- ↵
- P. Huang,
- X. Zhang,
- S. Ma,
- X. Huang

- ↵
- S. Dunatunga,
- K. Kamrin

- ↵
- J. Warren et al.

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Applied Physical Sciences