Ken Kamrin, David Henann MIT, Department of Mechanical Engineering Georg Koval

1. Constructing and verifying a three-dimensional nonlocal granular rheology Ken Kamrin, David Henann MIT, Department of Mechanical Engineering Georg Koval National Institute of Applied Sciences (INSA), Strasbourg

2. Stage 1: Local modeling. Static zones and flowing zones. Samadani & Kudrolli (2002) Kamrin, Rycroft, & Bazant (2007) (Experiment) (Discrete Element Method [DEM])

3. Inertial Flow Rheology: Simple Shear • Inertial number: • Inertial rheology: d = mean particle diameter s = the density of a grain packing fraction [Da Cruz et al 2005, Jop et al 2006] (Show Bagnold connection)

4. Dense Granular Continuum? (Blue crosshairs = Cauchy stress evecs) (Orange crosshairs = Strain-rate evecs) Rycroft, Kamrin, and Bazant (JMPS 2009) • Some evidence for continuum treatment: • Space-averaged fields vary smoothly when coarse-grained at width 5d • Certain deterministic relationships between fields appear in elements 5d wide

5. Developed packing fraction (5d scale) Tall silo Rycroft, Kamrin, and Bazant (JMPS 2009) Trap-door Wide silo

6. B Bt Local region To close the system: (1D picture) Deformed Reference Spring: Jiang-Liu granular elasticity law [Jiang and Liu, PRL, 2003] Dashpot: Inertial rheology [Jop et. al, Nature, 2006] Relaxed (3D version, Kroner-Lee decomposition) Constitutive Process Conservation of linear and angular momentum require

7. Mathematical Form Enforce: Frame indifference 2nd law of thermodynamics Coaxiality of flow and stress* “Mandel stress” Definitions: To close, must choose isotropic scalar functions  and dp with dp¸ 0.

8. Convert T(t) to M(t) using Fe(t)=F(t)Fp(t)-1. • Compute Dp(t) from M(t) using the flow rule. • Use Dp = Fp Fp-1 to solve for Fp(t+t). • Compute Fe(t+t) = F(t+t) Fp(t+t)-1. • Polar decompose Fe(t+t) and obtain Ee(t+t) = log Ue(t+t). • Compute M(t+t) from Ee(t+t) using the elasticity law. • Convert M(t+t) to T(t+t) using Fe(t+t). . Explicit Algorithm Create a user material subroutine Given input: F(t), Fp(t), T(t), and F(t+t) Compute output: Fp(t+t) and T(t+t)

9. (Y. Jiang and M. Liu, Phys. Rev. Lett., 2003) Granular Elasticity Law  FHertz /3/2 B is a stiffness modulus and  constrains the shear and compressive stress ratio. Under compression:

10. Results • Simulations performed using ABAQUS/Explicit finite element package. • Constitutive model applied as a user material subroutine (VUMAT). • All parameters in the model (6 in total) are from the papers of Jop et. al. and Jiang and Liu, except for d: K. Kamrin, Int. J Plasticity, (2010)

11. x y Rough Inclined Chute Velocity vector field: Experimental agreement: Elasto-plastic model Bingham fluid (no elasticity)P. Jop et. al. Nature (2006)

12. Wide Flat-Bottom Slab Silo Recall expt:

13. Cauchy stress field Eigendirections of stress: Theory Kamrin, Int. J Plasticity, (2010) DEM simulation Rycroft, Kamrin, and Bazant (JMPS 2009)

14. Flow prediction not “spreading” enough Theory Kamrin, Int. J Plasticity, (2010) Log of dp: DEM simulation Rycroft, Kamrin, and Bazant (JMPS 2009)

15. Flow prediction not “spreading” enough DEM data (Koval et al PRE 2009) (Width of shear band)/d Continuum model vwall Local model predictions are: Qualitatively ok: - Shear band at inner wall - Flow roughly invariant in the vertical direction Quantitatively flawed: - No sharp flow/no-flow interface. - Wall speed slows to zero, local law says shear-band width goes to zero. Continuum model Velocity field: x z G.D.R. Midi, Euro. Phys. Journ. E, (2004)

16. Stage 2: Fixing the problem.Accounting for size-effects • Past nonlocal granular flow approaches: • Self-activated process (Pouliquen & Forterre 2009) • Cosserat continuum (Mohan et.al 2002) • Ginzberg-Landau order parameter (Aronson & Tsimring 2001)

17. Deviating from a local flow law • Inertial number: • Local granular rheology: R/d R/d Pout 2R Koval et al. (PRE 2009) Vwall R Fix R/d, Vary Vwall I(x,y), (x,y) Vwall Fill with grains of diameter = d. d and Pout held fixed Fix Vwall , Vary R/d d = particle diameter s = the density of a grain [Da Cruz et al 2005, Jop et al 2006]

18. Nonlocal Fluidity Model Micro-level length-scale proportional to Existing theory for emulsions Extend to granular media Goyonet al., Nature (2007); Bocquetet al., PRL(2009) Define the “fluidity” (inverse viscosity): Define the “granular fluidity” : Local granular flow law: Local flow law (Herschel-Bulkley): Nonlocal add-on: Nonlocal law: Where bulk contribution vanishes, we can still flow even though we’re under the “yield criterion”.

19. Physical Picture Kinetic Elasto-Plastic (KEP) model Bocquet et al. (PRL 2009) Basic idea: Flow induces flow. Plastic dynamics are spatially cooperative. A local plastic rearrangement in box j causes an elastic redistribution of stress in material nearby. The stress field perturbation from the redistribution can induce material in box i nearby to undergo a plastic rearrangement. Granular Fluidity = g = “Relative susceptibility to flow”

20. Cooperativity Length Theoretical form for emulsions is Direct tests: From steady-flow DEM data in the 3 geometries (annular shear, vertical chute, shear w/ gravity): [Bocquet et al PRL 2009] Switching  Extract » from DEM using: Only new material constant is the nonlocal amplitude A. Local law constants all carry over. A=0.68 For our 2D DEM disks, we find:

21. Possible meaning of diverging length-scale M. van Hecke, Cond. Mat. (2009) Lois and Carlson, EPL, (2007) f>>fc f=fc+ Staron et al, PRL (2002) q=5o Pre-avalanche zone sizes for inclined plane flow: q=15o

22. Validating the nonlocal model Pout 2R Vwall R I(x,y), (x,y) Vary Vwall , Fix R/d R/d Vwall Check predictions against DEM data from Koval et al. (PRE 2009) in the annular couette cell. Kamrin & Koval PRL (2012) DEM (symbols) Nonlocal Rheology (lines) Local Rheology (- - -) Fix Vwall , Vary R/d R/d Wall slip? Perhaps part of joint BC.

23. Same model, new geometry: Vertical Chute G Pwall Velocity fields for 3 different gravity values: G DEM (symbols) Theory (lines) Kamrin & Koval PRL (2012)

24. Pwall Vwall Same model, new geometry: Shear with gravity y G H G Velocity fields for 3 different gravity values: DEM (symbols) Theory (lines) Kamrin & Koval PRL (2012)

25. Nonlocal rheology in 3D Solved simultaneously alongside Newton’s laws: Cooperativity length: Local rheology (standard inertial flow relation): 2 local material parameters: 2 grain parameters: 1 new nonlocal amplitude: For glass beads: + + Jop et al., JFM (2005) Our calibration for glass beads: A=0.48

26. Flow in the split-bottom Couette cell No existing continuum model has been able to rationalize flow fields in this geometry. Schematic: Filled with grains: grains H Inner section (fixed) Outer section rotates split van Hecke (2003, 2004, 2006)

27. Flow in the split-bottom Couette cell – surface flow The normalized revolution-rate: Viewed from the top down: van Hecke (2003, 2004, 2006)

28. Flow in the split-bottom Couette cell – surface flow Exp data from: Fenistein et al., Nature (2003) Flow fields on the top surface for five values of H: Rescaled: All shallow flow data collapses onto a near-perfect error function! Normalized radial coordinate

29. Nonlocal model in the split-bottom Couette cell • Custom-wrote an FEM subroutine to simulate the nonlocal model via ABAQUS User-Element (UEL). • Performed many sims varying H and d. Flow field on top surface Exp data from: Fenistein et al., Nature (2003) Simulated flow field in the plane: H=30mm H=10mm

30. Nonlocal model in the split-bottom Couette cell Define: Rc = Location of shear-band center w = Width of shear-band  = (r-Rc)/w = Normalized position Surface flow results for all 22 combinations of H and d tried: Collapse to error-function flow field Non-diffusive scaling of w with H

31. Nonlocal model in the split-bottom Couette cell - beneath the surface Exp data from: Fenistein et al., PRL (2004) Sub-surface shear-band width Sub-surface shear-band center The model describes sub-surface flow as well as surface flow.

32. Nonlocal model in the split-bottom Couette cell - Torque Normalized torque

33. Remove the inner wall and let H increase Rotation of the top-center point Exp data from: Fenistein et al., PRL (2006) The nonlocal model correctly captures the transition from shallow to deep behavior.

34. Existing 3D wall-shear data using glass beads Annular shear flow: Linear shear flow with gravity: Exp data: Losert et al., PRL (2000) Exp data: Siavoshi et al., PRE (2006) Same exact model and parameters used in split-bottom cases also captures these flows.

35. More on the cooperativity length Split-bottom cell Suppose we try different power laws in the cooperativity length formula: Annular shear Shear with gravity

36. Summary • Synthesized a local elasto-viscoplastic model by uniting the inertial flow rheology with a granular elastic response, and observed certain pros and cons. • Have proposed a nonlocal extension to the inertial flow law, by extending the theory of nonlocal fluidity to dry granular materials. Introduces only one new experimentally fit material constant. • Nonlocal model is demonstrated to reproduce the nontrivial size effect and rate-independence seen in slow flowing granular media. Vastly improves flow field prediction. • Have extended the model to 3D and demonstrated a strong agreement with experiments in many different flow geometries, including the split-bottom cell. Relevant papers: 2D prototype model: K. Kamrin and G. Koval, PRL (2012) 3D constitutive relation: D. Henann and K. Kamrin, PNAS (2013)

37. Important to do’s on dry flow modeling “Smaller is stronger”size effect Need more quantitative insights on fluidity BC’s (though much clarified by virtual power argument). BC’s likely responsible for observed nonlocal “thin-body” effects: Experiments (Silbert 2001) Nonlocal model w/ g=0 lower BC

38. Important to do’s on dry flow modeling “Smaller is stronger”size effect Focus on silo flow. Can we use our model to predict Beverloo constants that govern silo flow outflow rates? A famous size-effect problem. Beverloo constants Beverloo Correlation for drainage flow: Outflow rate Orifice size Particle size Q D D D=kd Q

39. Important to do’s on dry flow modeling Merge a transient model into the steady flow response: Easiest route: Let this term have transients per a critical-state-like theory Propose to try: (Anand-Gu form + Rate sensitivity) Let steady behavior approach inertial flow law [Rothenburg and Kruut IJSS (2004)]

40. Important to do’s on dry flow modeling Non-codirectionality of stress and strain-rate tensors Recall our usage of codirectionality (i.e. stress-deviator proportional to strain-rate tensor) However, data shows codirectionality is inexact in shear flows for two reasons: Slight non-coaxiality: Normal stress difference N2 (for 3D): (Shear) Depken et al., EPL (2007) (Black) Strain-rate eigen-directions (Blue) Stress eigen-directions

41. Important to do’s on dry flow modeling Material Point Method (MPM): Meshless method for continuum simulation. Useful for very large deformation plasticity problems. “Phantom” background mesh and a set of material point markers Our own preliminary MPM implementation: Great promise for full flow simulation [Z. Wieckowski. (2004)]