A Constitutive Model for Dense Granular Flows Based on

Article pubs.acs.org/IECR

A Constitutive Model for Dense Granular Flows Based on Microstructural Descriptors V. Vidyapati# and S. Subramaniam* Department of Mechanical Engineering, CoMFRE: Multiphase Flow Research & Education, Iowa State University, Ames, Iowa 50011, United States ABSTRACT: A constitutive model is developed to capture the complex rheological behavior of dense granular flows (solid volume fraction ranging from 0.45 to 0.62) in the quasi-static, intermediate, and inertial regimes. The principal contribution of this work is the development of a contact stress model (CSM) that is a statistical closure for the average contact stress experienced by particles, which is derived from a micromechanical model for the stress. This modeling approach naturally gives rise to the dependence of average contact stress on the average contact force and relevant descriptors of microstructure, which are the average coordination number and the fabric tensor. An expression for the average contact force is obtained from the contact force probability density function that has the same form in many granular flows. Appropriate closures for the coordination number and the fabric tensor are obtained by solving their respective modeled evolution equations proposed by Sun and Sundaresan [J. Fluid Mech. 2011, 682, 590−616]. CSM’s predictive capability is tested in homogeneous shear flow using discrete element model (DEM) simulation data corresponding to the quasi-static, intermediate, and inertial regimes. It is found that CSM’s predictions match DEM data very closely in all three regimes. Experiments performed by Tardos et al.10 reveal the existence of a third intermediate (transitional) regime that is characterized by σ ∝ γ̇n, where 0 < n < 2. These experiments also indicate that the intermediate regime is broad enough in the parameter space of solid volume fraction, particle friction coefficient, and the shear rate to require a separate model to capture its constitutive behavior. However, unlike the inertial (rapid flow) and the quasi-static regimes, the intermediate (transitional) regime still lacks a predictive constitutive model and has motivated many studies over the past decade.11−14 DEM simulations of canonical homogeneous shear flow14−16 as well as silo discharge2 confirm the existence of this intermediate regime of granular rheology. In a recent DEM study2 it is shown that all three regimes of granular rheology coexist in the silo discharge problem, and that the error in solid stress predicted by existing constitutive models correlates well with the spatial extent of the intermediate regime. Hence, a constitutive model that accurately captures the rheology of granular flow in all three regimes is necessary for predictive continuum simulations of granular flow. Figure 1 shows a schematic of a granular regime map with the corresponding constitutive behavior in each regime, where σ is the shear stress, γ̇ is the shear rate, and n is an exponent that relates the shear stress with the shear rate. Figure 1 also

1. INTRODUCTION Device-scale simulations of granular flow are required to accurately predict global flow characteristics in engineering applications such as hopper/silo discharge, chute flow, and dense-phase pneumatic conveying.1 In these large-scale simulations it is advantageous to use a continuum model that represents the average behavior of particles because the cost of simulating individual particle dynamics is prohibitively expensive, and often unnecessary for engineering purposes. These continuum simulations require a constitutive model that describes granular rheology. However, even for a simple flow such as silo discharge, continuum simulations using existing constitutive models predict discharge rates that differ considerably from three-dimensional (3D) DEM (discrete element method) simulations2 that compute individual particle dynamics. The DEM simulations are in good agreement with a correlation of experimental data,3 indicating that better constitutive models are needed in order for continuum models to predict global granular flow characteristics accurately. The rheology of dense granular flow is complex, and at least three different regimes have been identified. The inertial regime corresponds to rapid flow at relatively low solid volume fraction, and the kinetic theory for granular flow4−6 predicts a constitutive behavior in which the characteristic scale of stress increases as the square of the strain rate (σ ∝ γ̇2). At the other extreme, in the slow quasi-static flow regime, plasticity models applied to soil mechanics7,8 result in a stress that is independent of the applied shear rate (σ ≠ f(γ̇)). Both the inertial and quasistatic flow behavior are confirmed by DEM simulations.9 © 2016 American Chemical Society

Received: Revised: Accepted: Published: 10178

March 25, 2016 July 19, 2016 August 22, 2016 August 22, 2016 DOI: 10.1021/acs.iecr.6b01171 Ind. Eng. Chem. Res. 2016, 55, 10178−10190

Article

Industrial & Engineering Chemistry Research

stress−strain constitutive relation with density-dependent viscosity is proposed by Losert et al.21 In this model also the viscosity diverges when the density approaches the random close packing density of grains. Jop et al.11 proposed a constitutive relation inspired by the analogy between granular flows and viscoplastic fluids such as Bingham fluids for open surface flows (such as chute flow). In their work the effective friction coefficient, which is the ratio of shear stress to normal stress, is given as a function of the inertial number, I = γ̇d/(P/ρs)0.5, where d is the particle diameter, P is the pressure, and ρs is the particle density. Qualitative comparison of results for silo discharge using continuum simulations of the μ(I) − ν(I) model have been reported by Staron et al.22 using a variant of the original incompressible model in flows where dilatancy is important.23 The μ(I) − ν(I) model relies on the well-behaved limit of rigid contacts. Only in the limit of rigid contacts, in which deflections are so small that they are irrelevant, can one identify the inertial number I as the sole state parameter. This approach cannot be applied to situations in which grain elasticity (finite contact stiffness) matters, although Chialvo et al.16 have shown that departure from inertial number scalings is a direct result of particle softness, with a dimensionless shear rate characterizing the transition. The existence of the intermediate regime in which stresses vary with a different exponent as functions of the shear rate γ̇ is evidence for elastic deflections influencing the rheology. Therefore, the model of Jop et al.11 is not considered in this study. Experiments have shown that the average solid volume fraction is insufficient to describe the correct rheology of granular flows, and microstructural descriptors such as the coordination number and the fabric tensor have been measured using photoelastic particles.24,25 Also experiments in a 2D granular shear cell (GSC)26,27 as well as DEM simulations28 reveal that grain contacts in the intermediate regime are characterized by a mix of enduring solidlike and fluidlike contacts. In particular, these grain interactions are not determined by the solid volume fraction alone, but also depend on particle properties (such as particle friction coefficient and inelasticity) as well as on flow properties (such as the shear rate). Consequently, simple additive models where the total granular stress is decomposed into kinetic and frictional contributions are not able to capture the complex constitutive behavior in the intermediate (transitional) regime. As noted earlier, a recent study by the authors2 shows that the error incurred by existing constitutive models in predicting the solid stress for a silo discharge problem directly correlates with the spatial extent of the intermediate regime. Since most constitutive models in use are phenomenological, this observation motivates the development of a constitutive model for the intermediate regime that reflects the phase transition based on microstructural physical interactions between the grains. Vidyapati and Subramaniam14 showed that the refined order parameter (ROP) model that is based on an additive decomposition of the total granular stress into solidlike and fluidlike contributions, with each contribution being weighted by a function of the order parameter, is also incapable of reproducing intermediate regime rheology. By decomposing the total granular stress into contact and streaming contributions, Vidyapati and Subramaniam14 showed that the contact (virial) contribution to the stress dominates (>95%) in the intermediate regime, and follows the same scaling with strain rate as the total granular stress. However, the streaming contribution of the total stress always follows the

Figure 1. Schematic of granular regime map and their corresponding constitutive behavior in each regime.

shows the connection between the state/phase of granular matter and the rheology that is observed. Volfson and co-workers17 defined an order parameter for granular flows as the mesoscale space-time average of enduring solidlike contacts to the total contacts. They showed that the order parameter is null in the fluidlike inertial regime that is dominated by binary, instantaneous collisional contacts (adequately described by the kinetic theory of granular gases), whereas it takes the value unity in the solidlike quasistatic regime where enduring frictional contacts predominate. Using the same definition of the order parameter, Vidyapati and Subramaniam14 found that there is also a third stable phase with a liquidlike pair correlation in the intermediate regime. The development of an accurate constitutive model for dense granular flows remains an open and challenging problem partly due to the fact that the different constitutive behaviors of granular flow depend on both microscale properties (e.g., particle friction coefficient, coefficient of inelasticity) as well as macroscale properties (e.g., solid volume fraction and the shear rate). Of the three regimes, the intermediate regime where both collisional and enduring contacts between particles are important poses a significant challenge for constitutive model development. 14 The difficulty in developing accurate constitutive models of a postulated form is that it is difficult to know a priori which of the myriad particle and flow descriptors at different scales microscale, mesoscale, and macroscaleare needed to capture granular rheology. In the model developed in this work the relevant descriptors of microstructure emerge naturally from the derivation. Most constitutive models18−20 are based on an additive decomposition of the total granular stress as a weighted sum of kinetic and frictional contributions (σij = σijkin+σijfric), with the weight factor specified solely as a function of the solid volume fraction. The model for the kinetic contribution is readily available from some form of the kinetic theory and the model for the frictional contribution is postulated. A continuum theory for slow dense granular flows based on the so-called associated flow rule is proposed by Savage.18 This theory relates the shear stress and the strain rate in a plastic frictional system. Averaging strain-rate fluctuations yields a Bingham-like constitutive relation in which the shear stress has two contributions: a viscous part and a strain-rate independent part. According to this theory the stress and strain rate tensors are always coaxial. Furthermore, the theory also postulates that the viscosity diverges as the density approaches the close packing limit. A similar hydrodynamic model based on a Newtonian 10179

DOI: 10.1021/acs.iecr.6b01171 Ind. Eng. Chem. Res. 2016, 55, 10178−10190

Article

Industrial & Engineering Chemistry Research inertial scaling even in the intermediate regime. This result indicates that viscosity-based models may not be able to capture the correct granular rheology in the intermediate regime. This is because the shear viscosity is connected to the streaming part of the stress tensor in the kinetic theory, and it contributes less than 5% to the total granular stress in the intermediate and dense regimes. The ROP model (in both its KT (kinetic theory) and FSM (frictional stress model) variants14) fails because the solidlike stress contribution lacks the capability to reproduce the strain rate dependence found in intermediate regime rheology. These findings motivate the need for an accurate contact stress model in the intermediate regime of granular flow that properly accounts for the contact stress contribution to the total granular stress. In this work we develop a contact stress model (CSM) that is based on a statistical closure for the average contact stress experienced by a particle. This average contact stress is shown to depend on the average normal contact force that is obtained from the probability density function of normal contact force.29 The average contact stress also depends on microstructural descriptors such as the coordination number NCN and the fabric tensor R. Appropriate closures for the coordination number and the fabric tensor are obtained by solving their respective modeled evolution equations.30 To predict the total granular stress in the inertial regime, we couple the contact stress model (CSM) to the refined order parameter (ROP) model14 using the order parameter (OP) concept.17 Finally, the predictive capability of CSM is demonstrated by comparison with DEM data of homogeneous shear flow in different regimes of granular rheology: quasi-static, intermediate, and inertial. The following section presents the development of the CSM for dense granular flows. Appropriate closures that are needed for unclosed terms that appear in the CSM are presented in section 3. Section 4 details the DEM simulations of homogeneously sheared granular flow. Verification of closure terms are detailed in section 5. In section 6, we assess the performance of the CSM with DEM data for homogeneous shear flow in different regimes. Section 7 presents a comparative assessment of different constitutive models in the intermediate regime. The validity and consequences of the principal result are discussed in section 8 and conclusions are drawn in section 9.

c σαβ

1 [∑ m(i)vα′ (i)v′β(i)] V i

j,j≠i

(3)

Figure 2. Schematic of two particles i and j in contact with normal overlap δ(i)(j) and position vectors r(i), and r(j).

suffix notation is used with summation implied over repeated subscript indices, while particle indices are superscripted within parentheses and their summation is represented explicitly. The summation in the streaming stress (eq 2) runs over all the particles, while the summation in the contact stress (eq 3) runs over all the contacts in the averaging volume V. As noted earlier, recent DEM simulations by Vidyapati and Subramaniam14 reveal that for dense granular flows in the intermediate regime the contact contribution to the total granular stress dominates (contributing more than 95% of the total granular stress), and it has the same scaling with strain rate as the total granular stress. However, the streaming contribution always follows the inertial scaling (σ ∝ γ̇2) even in the intermediate regime. These results indicate that accurate modeling of the contact (virial) contribution of the granular stress is critical to capture the correct scaling of stress with strain rate in the intermediate regime. Note that eq 3 that is used to calculate the contact stress in DEM31,32 is essentially a stress estimator. A statistical contact stress model (CSM) based on this micromechanical stress estimator is obtained by taking the ensemble average of eq 3: CSM σαβ =

1 V

∑∑ i

j,j≠i

1 (i)(j) (i)(j) rα f β 2

1 ⟨Nc(V )⟩⟨rα(i)(j)f β(i)(j) ⟩ (4) 2V where ⟨Nc(V)⟩ is the average number of contacts in the sampling volume V, fβ(i)(j) is the contact force between particles i and j, and the ensemble average ⟨·⟩ denotes averaging over all multiparticle contacts. To completely specify the CSM model for average contact stress, the following quantities on the righthand side of eq 4 need to be modeled: (i) The average number of contacts ⟨Nc(V)⟩: This can be written in terms of the average (bulk) coordination number as follows:33 ≈

(1)

called the streaming and contact contributions, respectively. For soft spheres31,32 the expressions for the streaming and contact contributions to the stress tensors are given as follows. The streaming contribution arises from the flux of momentum associated with fluctuating particle velocity, and is identified as str σαβ =

∑

⎤ 1 (i)(j) (i)(j) ⎥ rα f β ⎥⎦ 2

where r(i)(j) is the separation vector pointing from the center of α particle j to the center of particle i, and fβ(i)(j) is the contact force exerted on particle i by particle j (see Figure 2). Standard tensor

2. CONTACT STRESS MODEL FOR DENSE GRANULAR FLOWS The stress tensor in granular media can be decomposed in two parts:31 σ = σ str + σ c

⎡ 1⎢ = ∑ V ⎢⎣ i

(2)

NCN ⟨N (V )⟩ (5) 2 where NCN is the average or bulk coordination number and ⟨N(V)⟩ is the average or expected number of particles in volume V. The average coordination number ⟨Nc(V )⟩ =

where m is the mass of a particle, v′α is the fluctuating velocity, i is a particle index, and the averaging volume is V. The contact contribution arises from interparticle interactions and is given by 10180


Article


The tangential component of the contact force is modeled by assuming that all contacts satisfy the Coulomb criterion at slippage: f t( i )(j) = μ f (i)(j) . The contact force can now be n written as

NCN is the average number of contacts per particle, and in DEM simulations it is computed as NCN =

∑i Nc(i) N − N1

(6)

f (i)(j) = f n(i)(j) (n(i)(j) + μt(i)(j))

N(i) c

where is number of contacts with the ith particle, N is the total number of particles in volume V, and N1 is the number of floaters (particles with less than one contact in a zero gravity environment). The coordination number is well-defined for all systems with N > N1. Besides it is expected that for very dilute granular flows where N1 approaches N, the CSM model will not be the appropriate rheological model for the granular stress but KT models will be more accurate. fβ(i)(j)⟩, the average (ii) The second required quantity is ⟨r(i)(j) α virial stress associated with a contact. Radjai et al.34 expressed ⟨r(i)(j) fβ(i)(j)⟩, the average virial stress α associated with a contact, as the integral over the joint probability density of forces and branch vectors projected on the shear plane. They decompose the contact force on the shear plane into its normal and tangential components, and the unit branch vector is parametrized by its orientation θ. The joint probability density and the normal and tangential contact force components are approximated by their lowest-order Fourier expansions in θ, with each expansion containing an anisotropy parameter and privileged direction. Their approach leads to unique insights into the effective friction coefficient that is approximated by one-half the sum of the anisotropy parameters. They show that the contact anisotropy (corresponding to the expansion of the joint probability density) is the principal cause of increase of effective friction with inertial number, which collapses several numerical and experimental data on dense sheared granular flows with different boundary conditions. They find that the contact anisotropy varies oppositely with the coordination number Z. The normal force anisotropy reflects force chains, which are increasingly destabilized by particle inertia causing it to decrease. The friction force anisotropy reflects friction mobilization and is responsible for force transmission. Apart from these interesting physical insights, their work shows the importance of contact anisotropy and coordination number in determining the effective friction, establishing the key particle-based parameters that determine macroscale behavior. The model presented in this paper takes a more engineering approach toward developing a better rheological model for the contact stress by recognizing the importance of contact anisotropy and coordination number as represented by statistics of the fabric tensor and average coordination number. 2.1. Average Contact Stress and Normal Contact Force. The contact force f(i)(j) on the ith particle can be represented in a coordinate system defined by the contact with particle j. This coordinate system is defined by a unit tangent vector t(i)(j) that lies in the contact tangent plane, and the unit vector orthogonal to it, called the contact normal n(i)(j).35 For spherical particles in contact the contact normal is the unit separation vector that lies along the line joining the centers of contacting particles, so n(i)(j) = r(i)(j)/|r(i)(j)|. Expressed in contact coordinates, the contact force f(i)(j) on particle i is f (i)(j) = (f n(i)(j) n(i)(j) + f t(i)(j) t(i)(j))

(8)

Substituting this expression for the contact force on a particle into the expression for ⟨r(i)(j) fβ(i)(j)⟩ in eq 4, and noting that α (i)(j) (i)(j) (i)(j) rα = r nα , shows the dependence of the average virial stress per contact on the normal contact force: ⟨rα(i)(j)f β(i)(j) ⟩ = ⟨rα(i)(j)f n(i)(j) (nα(i)(j)nβ(i)(j) + μnα(i)(j)t β(i)(j))⟩

(9)

The magnitude of the separation vector can be expressed in terms of the particle diameter d and the normal overlap δ(i)(j) between particles i and j in soft-sphere DEM as r(i)(j) = d − δ(i)(j). We also expect the average normal contact force to be linearly proportional to the extent of normal overlap: f (i)(j) ∼ n δ(i)(j). With these substitutions and neglecting terms of O(δ2) for small normal overlaps, the average contact stress in CSM becomes CSM σαβ =

1 nNCNd⟨fn (nαnβ + μnαtβ)⟩ 2

(10)

where n is the number density. Assuming that the normal contact force f n is independent of the orientation of contact normals (similar to the assumption in Radjai et al.34 who neglect force-fabric correlations in their form of the joint probability density), the expression for the average contact stress in CSM becomes CSM σαβ =

1 nNCNd⟨fn ⟩⟨(nαnβ + μnαtβ)⟩ 2

(11)

2.2. Fabric Tensor. The normalized distribution of contacts in granular media is generally described by a three-dimensional second-order fabric tensor35,36 R αβ = ⟨nαnβ ⟩

(12)

While the fabric tensor is formed by the covariance of the tensor product of unit normal vectors corresponding to particle contacts, the orthogonal fabric tensor Tαβ is formed by the covariance of the tensor product of the unit normal vector and the unit tangential vector at contact: Tαβ = ⟨nαtβ⟩

(13)

Using these expressions for the fabric and orthogonal fabric tensor in eq 11 one can express the CSM average contact stress as CSM σαβ =

1 nNCNd⟨fn ⟩(R αβ + μTαβ) 2

(14)

The orthogonal fabric tensor Tαβ is computed in DEM by defining a unit vector t as follows t=

ft |ft |

(15)

We computed the orthogonal fabric tensor Tαβ using data from our DEM simulations of homogeneously sheared granular flow. It is found that for solid volume fraction greater than 0.57 (corresponding to the relatively dense regime of granular flow) the orthogonal fabric tensor Tαβ is equal to the fabric tensor Rαβ (the difference in shear and normal components is within 8%,

(7) 10181


Article

Industrial & Engineering Chemistry Research dNCN = d 2|D|[(NCN,crit − NCN) + 7.5(ν − νc)0.5 ] dt

in other words within numerical error). Assuming this equality, the final CSM expression for the average contact stress is CSM σαβ

1 = nNCNd⟨fn ⟩R αβ(1 + μ) 2

where d2 = 5.6 is a material parameter is the modulus of the strain rate tensor. and νc are functions of the interparticle friction coefficient μ as described in Sun and Sundaresan.30 3.2. Closure for the Fabric Tensor. To model the fabric tensor we solve its evolution equation proposed by Sun and Sundaresan.30 The fabric tensor R is a microstructural quantity that describes the anisotropy of the contact distribution in granular media.35,36,38 The components of the fabric tensor can be calculated on the basis of particle contact information:

(16)

or in nondimensional form it reads: CSM σαβ d

kn

=

1 d2 nNCN ⟨fn ⟩R αβ(1 + μ) 2 kn

(17)

The rate-dependence in the CSM model given by eq 17 comes through the dependence of the fabric tensor as well as the coordination number on the shear rate. Evolution equations for both the fabric tensor R (eq 23) and the coordination number NCN (eq 18) contain the shear rate. In subsequent sections we use eq 17 to compute the average contact stress for a homogeneously sheared granular assembly and compare it with DEM results and other models. However, we first require closure relations for the average coordination number NCN, the fabric tensor Rαβ, and the average normal contact force ⟨f n⟩ to complete our specification of the contact stress model in eq 17.

R ij =

(22)

(23)

In eq 23, Ŕ = Ṙ + R·W − W·R, where W is the spin tensor, 1 W = 2 (∇u − (∇u)T ), and Ṙ denotes the material time deriv1

ative of R. In eq 23 |D| = 2 (DijDji) denotes the modulus of the strain rate tensor. Model constants c1 = −0.52, c2 = −2.8, and c3 = 100 are taken from Sun and Sundaresan.30 3.3. Average Normal Contact Force. Mueth et al.29 reported measurements of the distribution of normal forces exerted by granular material under uniaxial compression onto the interior surfaces of a confining vessel. These experiments using three-dimensional, random packings of monodisperse glass beads showed that this distribution is nearly uniform for forces below the mean force and decays exponentially for forces greater than the mean. The shape of the distribution was shown to be robust and unaffected by changes in the system preparation history or in the boundary conditions. They proposed the following empirical functional form for the distribution

(18)

where d1, d2, and d3 are material parameters in this microstructure evolution equation. In eq 18, R is the fabric tensor, D is the strain rate tensor based on the average solid velocity, S is 1 the deviatoric strain rate tensor S = D − 3 tr(D)I , and ν is the solid volume fraction. For a steady simple shear flow eq 18 reduces to following form

2

)

P(fn ) = a(1 − be−fn )e−βfn

(24)

where a = 3.0, b = 0.75, and β = 1.5 provides an excellent fit over the whole force range measured and is also consistent with simulation data. We find from DEM simulations that are reported in the following section that this form of the distribution holds even for homogeneously sheared granular flow although the coefficients are slightly different. We use this distribution with the coefficient values a = 2.43, b = 0.71, and β = 1.52 that best fit the DEM data for homogeneous shear to model the average normal contact force, which is obtained as

(19)

The function f(ν) dictates how the coordination number varies with the solid volume fraction in steady simple shear. This function has the following form:30 f (ν) = NCN,crit + β1(ν − νc)β2

∑ ninj c∈V

Ŕ = c1S + c 2|D|R + c3(R: S)R

dNCN = d1(R: S − χ |S|) + d 2|D|(f (ν) − NCN) + d3tr(D) dt

dNCN = d 2|D|(f (ν) − NCN) dt

1 Nc

where Rij is a symmetric second rank tensor, Nc is the number of contacts, ni and nj are the unit vectors corresponding to the contact vector from particle center to point of contact. Structural anisotropy can be easily related to the shear (xz component) of the fabric tensor R for simple shear flows (mean flow in the x direction and the shear gradient in the z direction). Sun and Sundaresan30 postulated the following evolution equation for the fabric:

3. CLOSURE RELATIONS FOR THE CONTACT STRESS MODEL In this section we describe the closures required for the average coordination number, the fabric tensor, and the average normal contact force. 3.1. Closure for the Average Coordination Number. The average coordination number is the average number of contacts per particle in the contact network. The average (or bulk) coordination number can be obtained by solving its modeled evolution equation proposed by Sun and Sundaresan.30 The coordination number characterizes the connectivity of a granular assembly and it is shown to be an important parameter that captures the phase transition in granular rheology.14 Sun and Sundraresan30 postulated the following evolution equation for the coordination number NCN:

(

(21)

1 and |D| = 2 (DijDji) The parameters NCN,crit

30

(20)

where model constants β1 = 7.5 and β2 = 0.5. In eq 20, NCN,crit and νc are the critical coordination number and critical solid volume fraction, respectively. NCN,crit varies from 4 to 6 as the particle friction coefficient changes from infinity to zero in three dimensions.37 Using eqs 19 and 20, one can write the evolution equation for the average coordination number in steady simple shear flow as

⟨fn ⟩ =

∫ fn P(fn ) dfn

(25)

4. DEM SIMULATIONS OF HOMOGENEOUSLY SHEARED GRANULAR FLOW To generate benchmark data for model assessment and to verify the proposed closures for CSM model inputs, we performed 10182


Article


Figure 3. Evolution of the average coordination number NCN (a) for ν = 0.61, μ = 0.1, k*= kn/(ρsd3γ̇2) = 1.0 × 105, e = 0.7 and (b) for ν = 0.62, μ = 0.1, k*= kn/(ρs d3γ̇2) = 1.0 × 105, e = 0.7. The solid line is the solution to eq 21, whereas the symbols are the data obtained from DEM simulations.

change, which raises questions about the suitability of this algorithm to study homogeneous time−dependent flows. This difficulty can be greatly alleviated through the use of the SLLOD algorithm. The SLLOD algorithm originates from ideas in nonequilibrium statistical mechanics39 where nonequilibrium flows such as shear flow are simulated by applying a force to the entire system (as opposed to simply moving the boundaries of the system faster or slower, as done in LeesEdwards). Although we do not study time-dependent shear in this article, the SLLOD algorithm was applied to all the simulations performed in this study to be consistent with other work.

DEM simulations of monodisperse, noncohesive spheres of diameter d and mass m subjected to homogeneous shear (where the stress is independent of position) over a range of solid volume fractions ν, particle friction coefficient μ, and the shear rate γ̇. A soft-sphere model is used in which particles interact via contact laws and friction only on contact. Since the realistic modeling of particle deformation is complicated, a simplified contact force model based on a linear spring− dashpot combination is used in this work.32 Details of the computational model used in the discrete element simulations are given in Appendix A. These constant-volume DEM simulations of sheared granular flow are performed in a cubical domain of side length L = 14d. The effect of system size is examined by varying the box length from 7d to 20d. It was found that the stress asymptotes once the box length exceeds 10d, consistent with estimates reported by Campbell.9 For all the simulations reported, the mass and diameter of the particles are set to 1, so the density of the particles is 6/π. The value of normal spring stiffness kn is set to 2 × 105 (mg/d units), which captures the general behavior of intermediate to high kn systems.32 The value of the coefficient of restitution e is chosen to be 0.7. All these simulations are performed with zero gravity. The integration time step Δt for all the simulations is selected to be tc/50, where tc is the binary collision time. This time step is shown to be sufficiently small to ensure temporal convergence.32 Simulations are run to a nondimensional time of γ̇t = 500, which is long enough to attain a statistically stationary solution.9 After reaching steady state the quantities are time-averaged over a time window corresponding to 200γ̇−1. In particular, since statistics in homogeneous shear flow are both stationary and homogeneous, the average stress is estimated by time-averaging the volumeaveraged estimate σαβ(t) given by eqs 1−3 as follows: ⟨σαβ⟩ ≈ ⟨σαβ⟩V , T =

1 T

∫t

t0+ T 0

σαβ (t ) dt

5. VERIFICATION OF CLOSURES FOR TERMS IN THE CONTACT STRESS MODEL We use the DEM simulations to verify the closures proposed for average coordination number and the fabric tensor in section 3. We also compare the distribution of normal contact force obtained from these DEM simulations of homogeneously sheared granular flow to that proposed by Mueth et al.29 from experiments of static granular packing. 5.1. Average Coordination Number. We compute the average coordination number from DEM data neglecting particles with zero contacts (floaters) or one contact (rattlers) as they do not participate in the contact network, consistent with the practice of other researchers.40 We use a first-order explicit scheme to march eq 21 in time. Figure 3 panels a and b show the evolution of the coordination number for a simple homogeneous shear flow for a solid volume of 0.61 and 0.62, respectively. The solid line is the solution to eq 21 and the symbols are the data obtained from DEM simulation. These figures show that average coordination number grows in time and attains a steady value at time γ̇t ≈ 2. As Figure 3 shows, the model for evolution of the average coordination number given by eq 21 closely follows the DEM data and reaches the same steady value. These results verify that the postulated form of the average coordination number evolution (eq 21) captures the correct steady value of the coordination number observed in DEM simulations of simple homogeneous shear flow. 5.2. Fabric Tensor. We solve eq 23 for the evolution of the fabric tensor and compare the (xz) component with that obtained from DEM data of homogeneously sheared granular flow (the mean axial velocity ⟨Ux⟩ is sheared in the z coordinate direction such that ∂⟨Ux⟩/∂ z = γ̇). Figure 4 shows the

(26)

with T = 200γ̇−1 and t0 = 500γ̇−1. These homogeneous shear simulations are performed with periodic boundary conditions in all directions (x,y,z), and uniform shear is generated in the domain using the “SLLOD” algorithm.39 The SLLOD algorithm39 is used in conjunction with the Lees-Edwards boundary condition to generate simple shear flows. The shearing motion induced by the Lees-Edwards boundary condition takes time to develop. Therefore, the flow would not be homogeneous immediately after a shear rate 10183


Article


Figure 4. Evolution of the xz component of the fabric tensor R for ν = 0.61, μ = 0.1, k* = kn/(ρsd3γ̇2) = 1.0 × 105 and e = 0.7 for homogeneous shear flow. The solid line is the solution to eq 23, and symbols are the data obtained from DEM simulations.

Figure 5. Probability density function of the normal contact force from DEM data with the fit to eq 24.

packings; that is, it exhibits a peak (plateau) for force less than mean (f < 1) and exponential decay for large forces ( f > 1). To compare with this experimental data, we fit eq 24 to our computational data and find the function with a = 2.43, b = 0.71, and β = 1.52 agrees well with the force distribution in the force networks as shown in Figure 5. We have extracted this force distribution for different values of solid volume fraction ν, particle friction coefficient μ, and the shear rates γ̇, but none of these variations are found to significantly influence the shape of this force distribution and it remains robust. Mueth et al.29 also noticed that this force distribution is a robust property of these granular packings based on their experimental data. Having verified all the required closures for the contact stress model, one can use eq 17 to predict contact stresses in homogeneous granular flows.

evolution of the xz component of the fabric tensor for a simple homogeneous shear flow for a solid volume fraction of 0.61. The solid line is the solution to eq 23 and the symbols are the data obtained from DEM simulations. This result shows that the postulated form of the fabric tensor evolution (eq 23) is capable of capturing the qualitative behavior observed in the DEM data, and its prediction of the steady value of the fabric tensor is within 10% of DEM. The slight underprediction of the fabric tensor component in Figure 4, occurs because the model for fabric evolution (eq 23) does not account for dependence on the solid volume fraction and particle friction coefficient. Table 1 shows the steady values of the fabric tensor obtained by Table 1. Variation of Fabric Tensor R with Nondimensional Shear Rate k*= kn/(ρs d3γ̇2) for ν = 0.61, μ = 0.1, and e = 0.7 k* = kn/(ρs d3γ̇2)

Rxz

× × × ×

−0.056 −0.048 −0.041 −0.035

2.5 1.0 1.0 1.0

104 105 107 109

6. CONTACT STRESS MODEL ASSESSMENT In this section we assess CSM’s capability to quantitatively predict granular rheology in the intermediate, quasi-static, and inertial regimes. We use a simple homogeneous shear problem for this purpose, and compare the magnitude of the shear stress predicted using the contact stress model (for a range of solid volume fractions ν, particle friction coefficients μ, and the shear rates γ̇) with DEM data. As noted earlier, the contact stress is the dominant contributor to the total granular stress (more than 95%) in the intermediate and quasi-static regime of granular flow.14 Hence, the capability of CSM to predict the contact stress is assessed in the intermediate and quasi-static regimes. However, in the inertial regime the streaming contribution to the total granular stress can be significant (up to 15% for a solid volume fraction of 0.45 with a particle friction coefficient of 0.1). Therefore, prediction of the total granular stress is required in this regime. To accomplish this, the CSM is coupled to the ROP model as described in section 6.3. These tests assess CSM’s performance in predicting the correct power-law dependence of the shear stress on the shear rate in different regimes of granular rheology. 6.1. Intermediate Regime. Figure 6a shows a logarithmic plot of elastic scaling of the shear component of the contact stress as a function of the nondimensional stiffness k* = kn/ ρsd3γ̇2 for a solid volume fraction of 0.58 with a particle friction coefficient of 1.0. In this scaling, the inertial regime dependence of stress on shear rate (σ ∝ γ̇2) appears as a line with slope m = −1, whereas the quasi-static regime where stress is independent of shear rate appears as a horizontal line. Intermediate regime rheology where σ ∝ γ̇n with 0 < n < 2 is characterized by

solving eq 23 for different values of nondimensional shear rate k*= kn/(ρs d3γ̇2). DEM studies by Vidyapati41 indicate that the shear rate has the biggest influence on the fabric tensor, while the fabric is independent of interparticle friction as long as it is greater than 0.1. This independency of fabric tensor on interparticle friction coefficient (when μ > 0.1) is also noted by Sun and Sundaresan.30 Nevertheless, Vidyapati and Subramaniam14 have also shown that the fabric is not very sensitive to changes in the order parameter (and hence, to the solid volume fraction and particle friction coefficient) when compared to other microstructural descriptors such as the average coordination number for a homogeneously sheared granular assembly. 5.3. Distribution of Normal Contact Force. The closure for the average normal contact force was specified using the probability density function (PDF) of the normal force. The PDF of this random variable is constructed using normal force data obtained from our 3D DEM simulations of homogeneous shear flows. In Figure 5 we show the resulting force distribution P(f) (where f = fn / fn is the normalized force). Figure 5 shows that about 70% of the contacts carry higher magnitude of normal force than the average normal force. The PDF shows the generic features of the force distribution in granular 10184


Article


Figure 6. Shear component of contact stress as a function of the nondimensional spring stiffness (inverse shear rate squared) k* = kn/(ρsd3γ̇2) in the intermediate regime (a) for ν = 0.58, μ = 1.0, e = 0.7; (b) for ν = 0.61, μ = 0.1, e = 0.7; and (c) for ν = 0.62, μ = 0.1, e = 0.7.

Both figures reveal that the CSM is capable of capturing the correct scaling of stress with shear rate (σ ∝ γ̇n, where 0 < n < 2) in the intermediate regime. The maximum difference between CSM and DEM is 25% in the higher shear rate regions (lower k* values). For higher values of solid volume fraction and k* (lower shear rate), the model predictions seem to capture the correct quasi-static behavior (where stress remains independent of the shear rate) as seen in Figure 6c. From Figure 6a−c, it is evident that CSM is capable of capturing the complex rheological behavior of granular flow in the intermediate regime. CSM’s excellent performance is attributed to the fact that it incorporates microstructural descriptors (such as the coordination number and the fabric tensor) in a constitutive modeling framework to capture structure-dependent rheology. 6.2. Quasi-static Regime. To assess the performance of the contact stress model in the quasi-static regime, we selected two test cases with solid volume fractions of 0.60 and 0.62, and an interparticle friction coefficient of 1.0. As noted earlier, in this regime also the contribution of streaming stress is negligible (less than 3%). Hence we compare the contact stress predicted using CSM with DEM data. In this regime, the stress remains independent of the shear rate (σ ≠ f(γ̇)) as previously shown by Campbell.9 For the cases in Figure 7a,b the critical solid volume fraction value is νc = 0.57,30 so ν > νc. Figure 7 panels a and b are logarithmic plots of the elastic scaling of the shear component of contact stress as a function of the shear rate for a solid volume fraction of 0.60 and 0.62, respectively. The particle friction coefficient selected for both of these test cases is unity, which ensures that these test cases correspond to the quasi-static regime. The shear stress predictions

stress versus shear rate dependence that lies between the lines with slope 0 and −1. It is noteworthy that some cases in the intermediate regime have a solid volume fraction greater than the critical value corresponding to jamming ν > νc (ν = 0.58, μ = 1.0; νc = 0.57 and ν = 0.62, μ = 0.1; νc = 0.615), whereas one case (ν = 0.61, μ = 0.1; νc = 0.615) has ν < νc. As noted earlier, values of critical volume fraction νc are obtained from Sun and Sundaresan30 corresponding to a specific value of interparticle friction μ. The role of critical solid volume fraction νc is also recognized by Campbell9 through elastic deformation of grains in which he identified the present intermediate regime as elastic−inertial and the present quasi-static regime as a elastic−quasi-static regime. In Figure 6a the shear component of contact stress obtained from CSM is shown by open diamonds, whereas the filled squares are the data from DEM simulations. Figure 6a shows that the contact stress behavior obtained from both CSM and DEM follows the intermediate scaling (σ ∝ γ̇n, where 0 < n < 2) of stress with the shear rate. Furthermore, the shear stress predicted by CSM closely follows the data obtained from the DEM simulations. Figure 6 panels b and c are similar plots of the shear component of the contact stress with shear rate, but for a solid volume fraction of 0.61 and 0.62, respectively, with a particle friction coefficient of 0.1. These simulation parameters are selected to ensure that they also correspond to the intermediate regime of granular flow. For the case in 6b the critical solid volume fraction value is νc = 0.615,30 so ν < νc. For the cases in 6c the critical solid volume fraction value is νc = 0.615, so ν > νc. This reveals that in the intermediate regime the solid volume fraction ν can either be higher or lower than its critical value νc. 10185


Article


Figure 7. Shear component of contact stress as a function of the nondimensional spring stiffness (inverse shear rate squared) k* = kn/(ρsd3γ̇2) in the quasi-static regime (a) for ν = 0.60, μ = 1.0, e = 0.7 and (b) for ν = 0.62, μ = 1.0, e = 0.7.

Figure 8. Shear component of the total granular stress as a function of the nondimensional spring stiffness (inverse shear rate squared) k* = kn/ (ρsd3γ̇2) in the inertial regime: (a) ν = 0.45, μ = 0.5, e = 0.7, and (b) ν = 0.53, μ = 0.5, e = 0.7.

where σfij is the fluidlike stress, σsij is the solidlike stress, σ0 = σsii/ 3(1 − α) is the scale of the stress, and α and β are model coefficients. The model coefficients (α and β) of the ROP model are specified as following polynomial fits of the order parameter ρ

obtained using CSM show rate-independent quasi-static behavior that closely follows the DEM data for both solid volume fraction values that were simulated. These results show that CSM can be used to successfully predict granular rheology even in the quasi-static regime. 6.3. Inertial Regime. In the inertial regime the streaming contribution of the total granular stress cannot be neglected. For example, the streaming stress can contribute up to 15% of the total granular stress for a solid volume fraction of 0.45 with an interparticle friction coefficient of 0.1. Hence, accurate prediction of the total granular stress is required in the inertial regime. However, CSM can only predict the contact stress. To predict the total granular stress, we couple CSM to the refined order parameter (ROP) model developed in our earlier work.14 The ROP model relies on an additive decomposition of the total granular stress into solidlike and fluidlike contributions, the relative magnitude of which depends on the order parameter (OP)ρ. The OP obeys a Ginzburg−Landau equation14,28 and characterizes the transition of the granular flow from enduring solidlike contacts to instantaneous fluidlike contacts. The model equations for the ROP model are14 σijf

= σ0{αδij + βbij}

σijs = σ0{(1 − α)δij + (1 − β)bij}

α = a + bρ + cρ2 + dρ3

(29)

where a = 1.0, b = −1.23, c = −0.31, and d = 0.54, and β = A + Bρ + Cρ2 + Dρ3

(30)

where A = 1.0, B = −1.69, C = 0.76, and D = −0.07. CSM relies on a different decomposition of the total granular stress into contact and streaming parts, since these contributions reflect the dependence of stress on shear rate more accurately.14 The two concepts are easily combined by recognizing that the OP is defined as the ratio of enduring solidlike contacts to total contacts. Therefore, we propose that the OP should also reflect the ratio of the solidlike stress to the contact stress as σijs = ρσijc

(27)

(31)

where σcij is the contact stress and ρ is the order parameter. The total granular stress in the CSM−ROP model is then obtained as follows:14

(28) 10186


Article

Industrial & Engineering Chemistry Research σij = =

s ⎤ σ0 ⎡ σij ⎢ + δij(β + α)⎥ (1 − β) ⎣ σ0 ⎦

⎤ σ0 ⎢ + δij(β + α)⎥ (1 − β) ⎣ σ0 ⎦

where viscosity is a function of the density as follows,

η = (νmax − ν)−1.75

⎡ ρσijc

(34)

Figure 9 shows that the shear stress predicted using this model fails to capture the correct scaling of the shear stress with the shear rate in the intermediate regime. As noted earlier, this model fails because it is a viscosity-based model. The viscosity is connected to the streaming contribution of the stress in the kinetic theory. However, this streaming part contributes less than 5% to the total granular stress in the intermediate and dense regime. Hence this model is not able to capture the correct scaling of stress with shear rate in this regime of granular flow. It should be noted that this model21 was proposed based on experimental data obtained from shear flow of granular material in a Couette geometry. ROP−KT: This is the constitutive model proposed by Vidyapati and Subramaniam,14 in which the ROP (refined order parameter) model is coupled with the kinetic theory of granular flows (KTGF)4 for the fluidlike stress. As seen in Figure 9 this model also fails to capture the correct trends of shear stress with the shear rate in the intermediate regime. The reason for this failure is attributed to the fact that the ROP−KT model assumes that the fluidlike stress contribution follows the kinetic theory closure even in the intermediate regime of flow. However, it is found that this assumption does not hold in the intermediate regime.14 ROP−FSM: In this constitutive model the ROP model is coupled with the frictional stress model (FSM) developed by Srivastava and Sundaresan20 for the frictional part of the total granular stress. The FSM model is used to compute the solidlike stress contribution σsij and then the ROP model is solved to obtain the total granular stress as follows:

(32)

where σ0 = σsii/3(1 − α) is the scale of the stress, and α and β are model coefficients that are functions of the OP, which is obtained by solving the Ginzburg−Landau equation as described in Vidyapati and Subramaniam.14 Figure 8 panels a and b show a logarithmic plot of the elastic scaling of the shear component of the total granular stress as a function of shear rate for two cases in the inertial regime: one for a solid volume fraction of 0.45, and the other for 0.53. As noted previously, in this scaling the stress values in the inertial regime where σ ∝ γ̇2 lie on a line with slope −1 (in Figure 8a,b the slope of the line is denoted by m). The shear component of the total granular stress obtained from the CSM−ROP model (cf. eq 32 where σcij is obtained using CSM) is shown by the black diamonds, whereas the filled squares show the data from DEM simulations. Both CSM−ROP and DEM follow inertial scaling (σ ∝ γ̇2) of the stress with the applied shear rate. The total granular stress predicted using the CSM−ROP model closely follows the DEM data in both cases, and on this basis we conclude that the CSM−ROP model is able to capture inertial regime rheology as well. It is also interesting to note that for these cases in Figure 8a,b the critical solid volume fraction value is νc = 0.58,30 hence ν < νc.

7. COMPARATIVE EVALUATION OF DIFFERENT CONSTITUTIVE MODELS IN THE INTERMEDIATE REGIME We assess the performance of the CSM−ROP model with other constitutive models in the intermediate regime. Figure 9

σij =

s ⎤ σ0 ⎡ σij ⎢ + δij(β + α)⎥ (1 − β) ⎣ σ0 ⎦

(35)

σsii/3(1

where σ0 = − α) . One should note that eq 35 diverges as the OP goes to zero (its fluidlike limit), which reflects that the ROP−FSM model for the total stress does not contain any information about the fluidlike stress. This frictional stress model is based on the critical state theory of soil mechanics. At the critical state the granular assembly deforms without any volume change and the frictional contribution of the stress is given by σ fric =I− pc (ν)

2 sin ϕ

S S: S

(36)

where the form for pc(ν) (critical state pressure) is taken from Johnson and Jackson19

Figure 9. Shear component of the total granular stress plotted with the nondimensional spring stiffness (inverse shear rate squared) k* = kn/(ρsd3γ̇2). Simulation parameters are ν = 0.62, μ = 0.1, e = 0.7.

⎧ (ν − νmin)r ⎪F s if ν > νmin pc (ν) = ⎨ (νmax − ν) ⎪ if ν ≤ νmin ⎩0

shows the shear component of the total granular stress plotted with the nondimensional spring stiffness k* (inverse shear rate squared) for a solid volume fraction of 0.62 with an interparticle friction coefficient of 0.1. The different constitutive models that are compared are listed below: Losert (2000): Losert et al.21 developed a constitutive model with density-dependent viscosity. The shear stress in this model is given as σxy = ηγ ̇ (33)

(37)

where F, r, and s are constants, taken from Srivastava and Sundaresan.20 As shown in Figure 9, the ROP−FSM model predicts stresses that are independent of the shear rate (a behavior characteristic of the quasi-static regime). However, the data obtained from the DEM simulations show a dependency of the shear stress on shear rate that is characteristic of the intermediate regime. The ROP−FSM results show that it is not simply a matter of modeling the fluidlike or solidlike parts of 10187


Article


volume fraction values that lie in the realm of validity of KT it might be better to switch to the ROP−KT model since it is expected to be more accurate in that regime.

the total granular stress. Rather, what is lacking is a fundamental description of the dependence of stress on the shear rate in the intermediate regime. CSM−ROP: CSM is the constitutive model for the contact stress developed in the present work, in which stress is linked to the microstructural descriptors such as the coordination number and the fabric tensor, which is key to capturing the structure-dependent rheology correctly. To obtain the total granular stress, CSM is coupled to the ROP model as described in section 6.3. As seen in Figure 9 the CSM−ROP model is able to capture the correct scaling of the shear stress with the shear rate even in the intermediate regime of flow. This comparison shows that only the CSM−ROP model is able to correctly predict the scaling of shear stress with shear rate in the intermediate regime of granular flow.

9. CONCLUSIONS A constitutive model is developed to capture the complex rheological behavior of dense granular flows in quasi-static, intermediate, and inertial regimes. This constitutive model is based on a new model for the contact stress in dense granular flow. This statistical closure naturally incorporates the relevant microstructure descriptors, which are the average coordination number and the fabric tensor. In granular rheology these descriptors characterize the evolution of microstructure, which is eventually responsible for the complex rheological behavior observed in different regimes. To complete the specification of the contact stress model (CSM), we solve modeled evolution equations for the coordination number and the fabric tensor that have been proposed by Sun and Sundaresan.30 The steady state solution of these evolution equations are verified using DEM data of homogeneous simple shear flow. The closure for mean normal contact force is obtained by constructing the probability density function (PDF) of the normal force using data from DEM simulations. In dense granular flows the universal nature of this contact force PDF allows simplification, because the force distribution remains almost independent of particle (e.g., friction coefficient) and flow (e.g., shear rate) properties. The predictive capability of the contact stress model is verified by comparing its predictions for shear stress in homogeneously sheared flow with DEM data for a range of solid volume fractions and friction coefficients in the intermediate and quasi-static regimes. It is shown that the contact stress model is capable of capturing the complex rheological behavior of granular flows in the intermediate and quasi-static regimes. It is also shown that this model can be successfully extended to the inertial regime (where the streaming contribution to the total granular stress can no longer be neglected) by coupling it with the ROP model14 through the order parameter concept.17 The reason for the better predictive capability of the contact stress model is attributed to the fact that it is derived as the statistical form of the micromechanical expression for stress, unlike other models that postulate the form of the constitutive law. Consequently, this model naturally includes the relevant microstructure descriptors (the average coordination number and the fabric tensor) in the constitutive model, which eventually determine the rheological behavior of granular flows.

8. DISCUSSION This model for the contact stress in granular flow has been developed on the basis of several simplifying assumptions. It is worthwhile to examine their validity and explore how extensions may be developed to relax them. Particles are assumed to be spherical in the theoretical development but particles are often nonspherical in practical applications. While the general theoretical development will still hold for arbitrary particle shape, the specific forms of the contact force PDF and the closure models for the average coordination number, fabric tensor, and order parameter will change. The assumption that the PDF of normal contact force is independent of the orientation tensors will probably not hold for nonspherical particles, and data from experiment or DEM will be needed to extend that aspect of the model. The assumption that the orthogonal fabric tensor equals the fabric tensor is again something that will need to be reassessed for nonspherical particles. The assumption of local statistical homogeneity on the scale of the particle diameter is necessary for this statistical model, but macroscale inhomogeneities can be accommodated by allowing the average coordination number, fabric tensor, and order parameter to vary in space and time. Appropriate boundary conditions at walls and inflow/outflow boundaries will need to be developed for these quantities in order for CSM to be extended to inhomogeneous flows. All contacts are probably not at slippage and the contacts may not all obey the Coulomb criterion, but more sophisticated extensions can be developed on the basis of DEM data, or if data from experiment becomes available. The reason CSM predicts the contact stress so well is partly because structure and rheology are not very tightly coupled in this problem. The microstructural descriptors of structure such as the average coordination number and the fabric tensor depend only on the shear rate, and they are actually decoupled in the simple shear problem that is considered in this work. Another factor that contributes to the model’s success is the universality of the contact force distribution in these flows. For these reasons CSM is remarkably accurate in capturing intermediate regime rheology, but it is noteworthy that the coefficients in the microstructure evolution and the PDF of normal contact force have constants that are specific to the simple homogeneous shear flow considered here. Nevertheless, CSM forms the basis for a new class of predictive constitutive models for granular rheology that are based on microstructural quantities. In the CSM−ROP model the solid stress is computed using CSM while the fluid stress is obtained through the ROP relations. For dilute granular flows that correspond to solid

■

APPENDIX A. DEM PARTICLE CONTACT MODEL For two contacting particles {i,j}, with radii {ai,aj} at positions {ri,rj}, with velocities {vi,vj} and angular velocities {ωi, ωj} (see Figure 2), the normal compression δij, relative normal velocity vnij, and relative tangential velocity vtij are32 δij = d0 − rij

(A.1)

vnij = (vij·nij)nij

(A.2)

vtij = vij − vnij − (aiωi + ajωj ) × nij

(A.3)

where d0 = ai + aj, rij = ri − rj, nij = rij/rij, with rij = |rij| and vij = vi − vj. Note that there is no sum over repeated indices. The rate of change of the elastic tangential displacement utij, set to zero at the initiation of a contact is, 10188


Article

Industrial & Engineering Chemistry Research du tij dt

= vtij −

The authors would also like to thank Pranav Shrotriya for a useful comment on the contact force decomposition.

(u tij ·vij)rij rij 2

■

(A.4)

The last term in eq A.4 arises from rigid body rotation around the contact point and ensures that utij always lies in the local tangent plane of contact. Normal and tangential forces acting on particle i are Fnij = f (δij/d)(knδij nij − γnmeff vnij)

(A.5)

Ftij = f (δij/d)( −kt u tij − γtmeff vtij)

(A.6)

where kn,t and γn,t are the spring stiffness and viscoelastic constants, respectively, and meff = mimj/(mi+mj) is the reduced mass of spheres with masses mi and mj. The corresponding contact force on particle j is simply given by Newton’s third law, that is, Fji = −Fij. The function f(δij/d) = 1 is for the linear spring−dashpot model, and f (δij/d) = δij/d is for Hertzian contacts with viscoelastic damping between spheres. Static friction is implemented by keeping track of the elastic shear displacement throughout the lifetime of a contact. The static yield criterion, characterized by a local particle friction coefficient μ, is modeled by truncating the magnitude of utij as necessary to satisfy |Ftij| < |μFnij|. Thus, the contact surfaces are treated as “sticking” when |Ftij| < |μ Fnij|, and as “slipping” when the yield criterion is satisfied. The amount of energy lost in collisions is characterized by the value of the coefficient of restitution, which is defined as the negative ratio of the particle normal velocity after collision to the velocity before collision. For the linear spring−dashpot model, the coefficient of normal restitution and contact time can be analytically obtained: en = exp( −γntc/2)

■

(A.7)

where the contact time tc is given by 2

−1/2

tc = π (k n /meff − γn /4)

(A.8)

The value of the spring constant should be large enough to avoid particle interpenetration, yet not so large as to require an unreasonably small simulation time step Δt, since an accurate simulation typically requires Δt ∼ tc/50.9 After the contact force is calculated, the equations of motion, which are ordinary differential equations, can be numerically integrated to get the particle trajectories.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

# V.V.: Corporate Engineering Technology Laboratory, The Procter & Gamble Company, 8256 Union Centre Boulevard., West Chester, Ohio 45069, United States.

■

Notes

The authors declare no competing financial interest.

■

NOMENCLATURE d = particle diameter e = particle restitution coefficient I = inertia number k* = nondimensional shear rate kn = particle normal stiffness coefficient kt = particle tangential stiffness coefficient m = particle mass N = total number of particles n = number density N1 = number of floaters NCN,crit = critical coordination number NCN = average coordination number Nc = number of contacts P = pressure Pc = critical state pressure t = time tc = binary collision time V = averaging volume D = strain rate tensor f = contact force g = acceleration due to gravity I = identity tensor n = unit normal vector R = fabric tensor r = separation vector S = deviatoric strain rate tensor t = unit tangent vector v′ = fluctuating velocity W = spin tensor Δt = time step for DEM simulations GREEK SYMBOLS δ = normal overlap γ̇ = shear rate η = viscosity μ = interparticle friction coefficient v = solid volume fraction vc = critical solid volume fraction vmax = maximum solid packing vmin = minimum frictional solid packing ρ = order parameter ρs = particle density σ = shear stress σc = contact contribution of stress tensor σfric = frictional contribution of stress tensor σf = fluidlike contribution of stress tensor σstr = streaming contribution of stress tensor σs = solidlike contribution of stress tensor σtot = total stress tensor REFERENCES

(1) Sundaresan, S. Some outstanding questions in handling of cohesionless particles. Powder Technol. 2001, 115, 2−7. (2) Vidyapati, V.; Subramaniam, S. Granular flow in silo discharge: Discrete element method simulations and model assessment. Ind. Eng. Chem. Res. 2013, 52, 13171−13182. (3) Beverloo, W. A.; Leniger, H. A.; van de Velde, J. The flow of granular solids through orifices. Chem. Eng. Sci. 1961, 15, 243−250.

ACKNOWLEDGMENTS The authors wish to acknowledge the financial support from U.S. department of Energy National Energy Technology Laboratory (Grant DE-FG26-07NT43070), administered under Advanced Coal Research at U.S. colleges and universities. 10189


Article

Industrial & Engineering Chemistry Research (4) Lun, C.; Savage, S.; Jeffrey, D.; Chepurniy, N. Kinetic theories for granular flow: inelastic particles in Couette flow and slightly inelastic particles in general flow fields. J. Fluid Mech. 1984, 140, 223−256. (5) Jenkins, J. T.; Savage, S. B. A theory for the rapid flow of identical, smooth, nearly elastic, spherical particles. J. Fluid Mech. 1983, 130, 187−202. (6) Goldhirsch, I. Rapid granular flows. Annu. Rev. Fluid Mech. 2003, 35, 267−293. (7) Nedderman, R. M. Static and Kninematics of Granular Material, 2nd. ed.; Cambridge University Press, 1992. (8) Schaeffer, D. Instability in the evolution equations describing incompressible granular flow. Journal of Differential Equations 1987, 66, 19−50. (9) Campbell, C. Granular shear flows at the elastic limit. J. Fluid Mech. 2002, 465, 261−291. (10) Tardos, I.; McNamara, S.; Talu, I. Slow and intermediate flow of a frictional bulk powder in the couette geometry. Powder Technol. 2003, 131, 23−39. (11) Jop, P.; Forterre, Y.; Pouliquen, O. A constitutive law for dense granular flows. Nature 2006, 441 (8), 727−730. (12) MiDi, G. D. R. On dense granular flows. Eur. Phys. J. E: Soft Matter Biol. Phys. 2004, 14, 341−365. (13) Vidyapati, V.; Langroudi, M.; Sun, J.; Sundaresan, S.; Tardos, G.; Subramaniam, S. Experimental and computational studies of dense granular flow: Transition from quasi-static to intermediate regime in a Couette shear device. Powder Technol. 2012, 220, 7−14. (14) Vidyapati, V.; Subramaniam, S. Granular rheology and phase transition: DEM simulations and order-parameter based constitutive model. Chem. Eng. Sci. 2012, 72, 20−34. (15) Ji, S.; Shen, H. Internal parameters and regime map for soft polydispersed granular material. J. Rheol. 2008, 52, 87−103. (16) Chialvo, S.; Sun, J.; Sundaresan, S. Bridging the rheology of granular flows in three regimes. Phys. Rev. E 2012, 85, 1−8. (17) Volfson, D.; Tsimring, L.; Aranson, I. Order parameter description of stationary partially fluidized shear granular flows. Phys. Rev. Lett. 2003, 90, 1−4. (18) Savage, S. Analyses of slow high-concentration flows of granular materials. J. Fluid Mech. 1998, 377, 1−77. (19) Johnson, P.; Jackson, R. Frictional collisional constitutive relations for granular materials, with application to plane shearing. J. Fluid Mech. 1987, 176, 67−98. (20) Srivastava, A.; Sundaresan, S. Analysis of a frictional-kinetic model for gas-particle flow. Powder Technol. 2003, 129, 72−85. (21) Losert, W.; Bocquet, L.; Lubensky, T.; Gollub, J. Particle dynamics in sheared granular matter. Phys. Rev. Lett. 2000, 85 (7), 1428−1431. (22) Staron, L.; Lagree, P. Y.; Popinet, S. Continuum simulation of the discharge of the granular silo. Eur. Phys. J. E: Soft Matter Biol. Phys. 2014, 37, 1−12. (23) da Cruz, F.; Emam, S.; Prochnow, M.; Roux, J. N.; Chevoir, F. Rheophysics of dense granular materials: Discrete simulation of plane shear flows. Phys. Rev. E 2005, 72, 021309. (24) Oda, M.; Konishi, J.; Nemat-Nasser, S. Some experimentally based fundamental results on the mechanical behaviour of granular materials. Geotechnique 1980, 30, 479−495. (25) Subhash, G.; Nasser, S. N.; Mehrabadi, M. M.; Shodj, H. M. Experimental investigation of fabric-stress relations in granular material. Mech. Mater. 1991, 11, 87−106. (26) McCarthy, J.; Jasti, V.; Marinack, M.; Higgs, C. Quantitative validation of the discrete element methid using an annular shear cell. Powder Technol. 2010, 203, 70−77. (27) Jasti, V.; Higgs, C. Experimental study of granular flows in a rough annular shear cell. Phys. Rev. E 2008, 78, 041306. (28) Volfson, D.; Tsimring, L.; Aranson, I. Partially fluidized shear granular flows: Continuum theory and molecular dynamics simulations. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2003, 68, 021301.

(29) Mueth, D.; Jaeger, H.; Nagel, S. Force distribution in a granular medium. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 1998, 57 (3), 3164. (30) Sun, J.; Sundaresan, S. A constitutive model with microstructure evolution for flow of rate-independent granular materials. J. Fluid Mech. 2011, 682, 590−616. (31) Luding, S., Hinirichsen, H., Wolf, D. E. The Physics of Granular Media; Wiley-VCH Verlag: New York, 2004. (32) Silbert, L.; Ertas, D.; Grest, G.; Halsey, T.; Levine, D.; Plimpton, S. J. Granular flow down an inclined plane: Bagnold scaling and rheology. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 2001, 64, 051302. (33) Zhang, H.; Makse, H. Jamming transition in emulsions and granular materials. Phys. Rev. E 2005, 72, 011301. (34) Azema, E.; Radjai, F. Internal structure of inertial granular flows. Phys. Rev. Lett. 2014, 112, 078001. (35) Bathurst, R.; Rothenburg, L. Observation on stress force fabric relationships in idealized granular material. Mech. Mater. 1990, 9, 65− 80. (36) Cowin, S. Anisotropic poroelasticity: fabric tensor formulation. Mech. Mater. 2004, 36, 665−667. (37) Song, C.; Wang, P.; Makse, H. A. A phase diagram for jammed matter. Nature 2008, 453, 629−632. (38) Radjai, F.; Wolf, D.; Jean, M.; Moreau, J. J. Bimodal character of stress transmission in granular packings. Phys. Rev. Lett. 1998, 80 (1), 61−64. (39) Evans, D. J., Morriss, G. P. Statistical Mechanics of Nonequillibrium Liquids, 2nd. ed.; Academia Press, 1990. (40) Shundyak, K.; van Hecke, M.; van Saarloos, W. Force mobilization and generalized isostaticity in jammed packings of frictional grains. Phys. Rev. E 2007, 75, 010301. (41) Vidyapati, V. Constitutive modeling of dense granular flow based on discrete element method simulations. Ph.D. Thesis, Iowa State University, 2012.

10190


A Constitutive Model for Dense Granular Flows Based on

Recommend Documents